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Bridge  piers  are  often  subjected  to  lateral  loading  that  is  not  neglectible 

when  compared  to  vertical  loads.  Such  loading  conditions  may  include  wind,  water  and 

earthquakes.  In  order  to  develop  a  time  domain  analysis  for  the  nonlinear  dynamic 

response  of  piers   and  their  foundations  the  computer  program   Florida-Pier  was 

modified.    Florida-Pier    is    a   nonlinear   finite-element   program,    developed    at   the 

University  of  Florida,  designed  for  analyzing  bridge  pier  structures  composed  of  pier 

columns  and  cap  supported  on  a  pile  cap  and  piles  with  nonlinear  soil.  The  program  was 

developed  in  conjunction  with  the  Florida  Department  of  Transportation  (FDOT) 

structures  division.  The  piers  and  the  piles  are  modeled  as  nonlinear  3D  beam  discrete 

elements.  These  elements  use  the  true  material  stress-strain  curves  (steel  and  concrete) 

to  develop  its  behavior  and  stiffness  modeling.  Nonlinear  dynamic  capability  was  added 

to  these  elements  by  adding  to  the  stress-strain  curves  the  ability  to  represent  loading, 

unloading  and  reloading  behavior,  typical  of  dynamic  loading.  In  addition,  a  mass 


matrix  for  this  element  was  also  implemented.  The  mass  matrix  for  the  piles  cap  was 
also  developed.  The  cap  is  modeled  as  linear  shell  elements.  The  stiffness  matrix  for  this 
type  of  element  is  considered  to  remain  unchanged  during  the  dynamic  analysis.  A 
nonlinear  model  was  also  added  for  the  soil.  The  p-y  springs  generated  by  the  program 
now  have  the  capability  of  loading,  unloading  and  reloading,  just  like  the  steel  and 
concrete.  Such  models  can  then  be  applied  to  dynamic  analysis  of  reinforced  concrete 
structures  subjected  to  seismic,  impulsive,  or  wind  loads.  In  addition,  being  the  current 
state-of-practice,  modal  analysis  was  also  implemented.  Seismic  load  can  now  be 
applied  to  the  linear  structure  considering  the  nonlinear  behavior  of  the  foundation.  A 
description  of  the  analysis  used  to  model  the  nonlinear  dynamic  behavior  of  bridge 
piers,  as  well  as  the  implemented  nonlinear  material  behavior,  is  presented  in  this 
dissertation. 


CHAPTER  1 
INTRODUCTION 

Background 

Extensive  research  efforts  have  been  directed  to  the  nonlinear  response  of 
structures  subjected  to  extreme  load  events.  These  extreme  load  events  could  be  an 
earthquake  or  hurricane  for  a  building,  a  ship  impact  for  a  bridge,  or  the  effects  of  waves 
and  wind  action  for  offshore  oil  platforms.  Traditionally,  large  factors  of  safety  have 
been  used  in  such  cases,  resulting  in  over-conservative  design  and  cost  ineffectiveness. 
On  the  other  hand,  an  unsafe  design  could  result  in  catastrophic  human  and  economic 
losses.  Because  a  more  sophisticated  nonlinear  dynamic  analysis  is  computationally 
expensive,  these  structures  are  designed  using  factored  static  loads  to  account  for  the 
dynamic  effects.  This  procedure  is  acceptable  for  very  low  frequency  vibrations, 
however  the  introduction  of  non-linearity,  damping,  and  pile-soil  interaction  during 
transient  loading  may  significantly  alter  the  response. 

Because  in  recent  years  computers  have  become  much  faster  and  cheaper,  it  has 
become  possible  to  consider,  and  consequently  to  study,  the  dynamic  nonlinear  behavior 
of  structures  considering  many  factors  neglected  in  the  past.  In  this  dissertation,  the 
computer  program  Florida  Pier  (Hoit  et  al.,  1996),  which  will  be  referred  simply  as 
FLPIER  from  now  on,  has  been  modified  to  allow  the  nonlinear  dynamic  analysis  of 
bridge  piers.  FLPIER  is  a  computer  program  based  on  the  Finite  Element  Method 


developed  by  Drs.  Hoit,  Mcvay,  and  Hays  at  the  University  of  Florida  for  the  nonlinear 
static  analysis  of  bridge  piers.  Nonlinear  aspects  of  structural  analysis,  such  as  material 
and  geometric  non-linearity,  as  well  as  structure-soil  interaction  can  be  incorporated  into 
the  analysis  leading  to  more  accurate  results.  It  can  model  all  the  components  of  a 
bridge  pier  and  its  foundation,  such  as  pier,  pier  cap,  piles  cap,  piles,  and  soil,  as  shown 
in  Fig. 1-1.  The  pier,  pier  cap,  and  piles  can  be  represented  using  discrete  elements  that 
can  incorporate  the  effects  of  material  and  geometric  nonlinear  behavior.  More  details 
about  the  discrete  element  are  found  in  Chapter  3.  The  piles  cap  is  modeled  as  linear  9- 
node  shell  elements.  The  lateral  soil  resistance  is  modeled  as  nonlinear  p-y  springs, 
while  the  axial  resistance  is  modeled  as  nonlinear  t-z  springs. 


Bridge  Springs 


Substructure 


Piles 


Soil  layers 


Fig.  1  - 1 .  Bridge  pier  components 


FLPIER  is  now  used  by  many  DOTs  throughout  the  United  States  because  of  its 
reliability  and  ease  of  use.  Unlike  other  general  Finite  Element  programs,  like  ADINA 
and  SAP,  where  modeling  and  analyzing  can  be  time  consuming,  in  FLPIER  it  is  easy 
and  fast  for  the  user  to  perform  these  tasks  thanks  to  a  user  friendly  interface  for  model 
generation.  The  modification  of  soil  or  structure  parameters  in  the  model  is  not  difficult 
either.  The  results  can  also  be  seen  through  a  graphic  interface  that  is  currently  being 
updated.  In  the  modified  dynamics  version,  resulting  from  this  research,  speed  and  ease 
were  maintained,  allowing  the  user  to  easily  perform  the  nonlinear  dynamic  analysis  and 
modify  parameters  in  the  soil  or  structure  if  necessary.  Although  the  program  is  more 
suitable  for  the  analysis  of  bridge  piers,  other  types  of  structures  can  also  be  modeled. 

The  new  contribution  for  the  field  is  the  proposed  concrete  model.  This  model 
was  implemented  in  the  FLPIER  code  to  allow  the  nonlinear  dynamic  analysis  of 
reinforced  concrete  sections. 

Literature  Review 

Over  the  last  years  different  analytical  models  have  been  proposed  for  the 

analysis  of  reinforced  concrete  structures.  Models  for  these  types  of  structures,  which 

are  under  primarily  flexural  and  axial  loads,  can  be  classified  as: 

(i)  Simple  or  lumped  models. 

(ii)  Discrete  models. 

(iii)  Fiber  models. 

(iv)  Finite  element  models. 

Single-degree-of-freedom  (SDOF)  models  belong  to  the  first  class  of  analytical 

models.  In  this  class  of  models  it  is  assumed  that  the  structure's  response  to  an 

earthquake  is  dominated  by  its  first  natural  frequency,  allowing  the  system  to  be 


4 

represented  as  a  SDOF  system  with  lumped  mass  and  stiffness  properties  (Crandall 
(1956),  Craig  (1981),  Paz  (1985),  and  Chopra  (1995)).  A  more  general  representation 
for  multi-degree-of  freedom  systems  (MDOF)  is  derived  using  the  concept  of  shear 
building.  In  this  model  the  stiffness  of  each  story  is  represented  by  nonlinear  springs, 
and  the  beams  are  considered  to  be  infinitely  rigid.  Despite  its  simplicity  and 
satisfactory  performance  in  predicting  the  maximum  response,  this  class  of  models  does 
not  provide  enough  data  for  more  detailed  seismic  analysis.  Furthermore  for  more 
complicated  frames  the  assumption  that  the  beams  are  infinitely  rigid  may  not  be 
correct. 

In  the  second  class  of  analytical  models,  the  discrete  models,  there  is  a 
correspondence  between  the  analytical  model  and  the  actual  structure.  In  such  models,  a 
linear  elastic  element  and  a  nonlinear  spring  represent  the  structural  elements.  The  most 
common  case  is  that  of  a  nonlinear  spring  attached  to  both  ends  of  a  linear  beam 
element.  Atalay  (1975),  Clough  (1966),  Nakata  et  al.  (1978),  Park  (1984),  and  Takeda 
(1970),  among  others,  have  extensively  used  this  class  of  models  to  analyze  the  behavior 
of  reinforced  concrete  structures.  In  these  models  a  set  of  predefined  rules  defines  the 
hysteretic  behavior  of  the  nonlinear  springs.  These  rules  are  usually  obtained  from 
laboratory  experiments  with  real  scale  specimens.  It  is  mainly  the  difference  between 
these  rules  that  distinguishes  the  models.  Although  these  models  give  satisfactory 
results,  its  main  disadvantage  is  the  fact  that  the  nonlinear  spring's  rules  are  based  on 
experiments  that  may  not  correspond  to  the  actual  structural  member,  or  type  of  loading, 
that  they  are  representing. 


The  fiber  models  have  been  used  in  the  study  of  reinforced  concrete  (Ala 
Saadeghvaziri  (1997),  Hajjar  et  al.  (1998),  Park  et  al.  (1972),  and  Zeris  and  Mahin 
(1991 -a),  Zeris  and  Mahin  (1991-b))  and  steel  members  (Baron  and  Venkatesan  (1969), 
Chen  and  Atsuta  (1973)).  These  models  are  based  on  the  finite  element  approach,  and 
are  better  suited  for  members  and  structures  under  complex  loading  histories.  In  these 
models  the  cross-section  is  divided  into  segments.  Each  segment  can  then  be  divided  in 
one  or  more  fibers.  Each  fiber  is  assumed  to  obey  a  uniaxial  stress-strain  relationship. 
From  the  integration  of  the  stresses  of  each  fiber  over  the  cross  section,  the  element 
forces  can  be  calculated,  and  from  the  evaluation  of  the  stiffness  of  each  fiber  the 
overall  element  stiffness  can  be  obtained.  Once  the  element  forces  and  stiffness  are 
obtained  the  analysis  is  carried  out  using  standard  Finite  Element  Method  procedures. 
Therefore  only  the  stress-strain  relationships  for  concrete  and  reinforcing  steel  in  the 
case  of  reinforced  concrete  sections,  or  steel,  in  the  case  of  steel  sections,  are  necessary 
to  describe  the  properties  of  each  section  of  the  element.  This  makes  these  models  very 
effective  under  complex  loads.  The  main  difference  among  all  the  fiber  models  are  the 
rules  adopted  for  the  uniaxial  behavior  of  the  different  materials  that  make  the  cross- 
section.  In  the  case  of  most  civil  engineering  structures,  steel  and  concrete,  but  other 
materials  can  also  be  used  if  the  stress-strain  relationships  are  known. 

In  the  specific  case  of  concrete  the  model's  backbone  is  the  envelope  curve 
obtained  from  a  monotonic  test.  This  curve  limits  the  concrete  stresses  in  any  loading 
phase.  In  some  models  the  compression  envelope  curve  for  concrete  is  represented  by 
the  well-known  Hognestad  parabola  (Ala  Saadeghvaziri  (1997),  Park  et  al.  (1972)). 
Another   approach   is   to   use   multilinear   curves   to   define   concrete    behavior   in 


compression  (Zeris  and  Mahin  (1991 -a),  Zeris  and  Mahin  (1991-b)).  In  the  case  of 
dynamic  analysis,  the  unloading  and  reloading  rules  are  particular  for  each  author.  In  the 
case  of  Ala  Saadeghvaziri  (1997),  and  Zeris  and  Mahin  (1991 -a),  Zeris  and  Mahin 
(1991-b),  unloading  under  compressive  stress  has  a  slope  equal  to  the  initial  Young's 
modulus  of  the  material.  However  Hajjar  et  al.  (1998)  and  Park  et  al.  (1972),  have  their 
own  more  complicated  expressions  for  the  unloading  curve.  The  reloading  curves  are 
very  specific  for  each  model,  and  the  reader  is  referred  to  the  mentioned  references  for 
more  information.  The  tension  strength  of  concrete  can  be  neglected  (Zeris  and  Mahin 
(1991 -a),  Zeris  and  Mahin  (1991-b)),  or  assumed  to  be  equal  to  the  concrete  tensile 
strength,  (Ala  Saadeghvaziri  (1997),  Park  et  al.  (1972))  in  which  case  its  slope  is 
assumed  to  be  equal  to  the  initial  slope  of  the  compression  side.  The  unloading  and 
reloading  criteria  are  again  specific  for  each  model,  and  the  reader  is  referred  to  the 
above  references  for  more  information.  Confinement,  strain-rate,  and  stiffness 
degradation  effects  are  also  particular  for  each  model. 

In  the  case  of  steel,  no  distinction  is  made  between  steel  sections  and 
reinforcement  steel.  The  rules  are  valid  for  both  cases.  For  the  stress-strain  curves,  the 
Baushinger's  effect  can  be  considered  (Park  et  al.  (1972),  Baron  and  Venkatesan 
(1969)),  or  ignored,  in  which  case  a  bilinear  or  tri-linear  relationship  (elastoplastic  with 
kinematic  or  isotropic  hardening)  is  used  (Ala  Saadeghvaziri  (1997),  Chen  and  Atsuta 
(1973),  Zeris  and  Mahin  (1991-a),  Zeris  and  Mahin  (1991-b)). 

The  major  disadvantage  of  these  models  is  that  they  are  computationally  very 
expensive,  but  with  the  recent  advances  in  computer  technology,  this  class  of  models 
has  become  more  popular  because  of  its  versatility.  This  is  the  analytical  model  used  in 


this  work.  FLPIER  already  incorporates  a  nonlinear  discrete  element,  which  uses  fiber 
modeling  at  two  points  along  the  element's  length  to  characterize  its  nonlinear  behavior. 
New  stress-strain  curves  were  introduced  to  allow  the  nonlinear  discrete  element  to 
perform  dynamic  analysis.  The  details  can  be  found  in  Chapter  3. 

The  last  class  of  analytical  methods  is  the  finite  element  method.  In  this  class  of 
methods  different  elements  are  used  to  represent  the  structural  members,  such  as  truss 
members  to  represent  the  reinforcing  steel,  and  plane  stress  elements  to  represent  the 
concrete.  The  cracking  typical  of  concrete  represents  a  computational  difficulty  for  these 
models,  requiring  the  development  of  more  sophisticated  elements.  The  development  of 
such  elements  is  a  challenge  based  on  complex  elasticity  and  plasticity  theories.  Like 
fiber  modeling,  this  class  of  methods  is  computationally  expensive  and  time  consuming. 

Limitations 

It  is  very  difficult  for  a  model  to  incorporate  all  the  aspects  inherent  to  nonlinear 
dynamic  analysis,  and  the  model  presented  here  is  no  exception.  The  first  limitation  is 
the  fact  that  all  the  theory  developed  in  Chapter  2  for  nonlinear  dynamic  analysis  is 
based  on  small  displacement  theory.  The  second  limitation  comes  from  the  fact  that  the 
effect  of  shear  deformation  is  not  included  in  the  constitutive  models.  It  was  also  not 
included  in  the  original  derivation  (Hoit  et  al.,  1996).  Problems  of  local  buckling  are 
also  outside  the  scope  of  this  work. 

Organization 

In  Chapter  2  all  the  theory  necessary  for  the  formulation  of  the  problem  is 
presented.  In  Chapter  3  the  derivation  of  the  discrete  element  is  described  and  the 
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adopted  constitutive  models  for  concrete  and  steel  are  presented.  Then  Chapter  4 
discusses  the  actual  state-of-design  procedure  for  modal  analysis.  In  Chapter  5  the 
concept  of  multiple  support  excitation  is  introduced.  In  Chapter  6  the  soil  structure 
interaction  and  the  dynamic  soil  behavior  are  explained.  In  Chapter  7  the  response 
predicted  by  FLPIER  is  compared  to  various  literature  results.  Finally  in  Chapter  8  the 
conclusions  and  suggestions  for  future  work  are  discussed. 


CHAPTER  2 
NONLINEAR  DYNAMIC  ANALYSIS 


Theory 


In  a  static  problem  the  frequency  of  the  excitation  applied  to  the  structure  is  less 
than  one  third  of  the  structure's  lowest  natural  frequency.  In  this  case  the  effects  of 
inertia  can  be  neglected  and  the  problem  is  called  quasistatic.  For  such  problems  the 
static  equations  [£]{£>}  =  {.R}are  sufficiently  accurate  to  model  the  response,  even 
though  the  loads  {R\ ,  and  displacements  {£>} ,  vary  (slowly)  with  time.  The  static  loads 
{ R]  may  result  from  surface  loads  and  /or  body  forces. 

On  the  other  hand  if  the  excitation  frequencies  are  higher  than  noted  above  or  if 
the  structure  vibrates  freely,  the  inertia  effects  must  be  considered  in  the  analysis.  The 
inertia  effects  are  accounted  for  by  the  mass  matrix,  written  as  [m]  for  an  element  and 
[M]  for  a  structure,  which  is  a  discrete  representation  of  the  continuous  distribution  of 
mass  in  a  structure.  The  effects  of  damping,  if  important,  are  accounted  for  by  the 
damping  matrices  [c]  and  [C]. 

The  dynamics  problems  can  be  categorized  as  either  wave  propagation  problems 
or  structural  dynamic  problems.  In  wave  propagation  problems  the  loading  is  often  an 
impact  or  an  explosive  blast.  The  excitation  and  the  structural  response  are  rich  in  high 
frequencies.  In  such  problems  we  are  usually  interested  in  the  effects  of  stress  waves. 
Thus  the  time  duration  of  analysis  is  usually  short  and  is  typically  of  the  order  of  a  wave 
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transversal  time  across  a  structure.  A  problem  that  is  not  a  wave  propagation  problem, 
but  for  which  inertia  is  important,  is  called  a  structural  dynamics  problem.  In  this 
category,  the  frequency  of  excitation  is  usually  of  the  same  order  as  the  structure's 
lowest  natural  frequencies  of  vibration. 

A  typical  example  of  a  wave  propagation  problem  would  be  that  of  analyzing  the 
stresses  in  a  pile  when  it  is  grounded.  It  is  important  not  to  exceed  the  allowable  stresses 
in  order  not  to  damage  the  pile.  Earthquake  analysis  of  structures  is  a  typical  structural 
dynamics  problem,  where  the  inertia  forces  govern  the  response  of  the  structure. 

Problems  of  structural  dynamics  can  still  be  subdivided  into  two  broad 
classifications.  In  the  first  one,  we  are  interested  in  the  natural  frequencies  of  vibration 
and  the  corresponding  mode  shapes.  Usually,  we  want  to  compare  natural  frequencies  of 
the  structure  with  frequencies  of  excitation.  In  design,  it  is  usually  desirable  to  assure 
that  these  frequencies  are  well  separated.  In  the  second  classification,  we  want  to  know 
how  a  structure  moves  with  time  under  prescribed  loads,  like  under  impacts,  blasts  or 
wind  loads,  and/or  motions  of  its  supports,  like  in  the  case  of  an  earthquake.  We  are 
interested  in  the  time-history  analysis.  The  two  most  popular  methods  of  dynamic 
analysis  are  modal  methods  and  direct  integration  methods. 

Equations  of  Motion.  Mass.  and  Damping  Matrices 

The  equations  that  govern  the  dynamic  response  of  a  structure  will  be  derived  by 
requiring  the  work  of  external  forces  to  be  absorbed  by  the  work  of  internal,  inertial  and 
viscous  forces,  for  any  small  kinematically  admissible  motion  (i.e.,  any  small  motion 
that  satisfies  both  compatibility  and  essential  boundary  conditions).  For  a  single 
element,  this  work  balance  becomes 


j{5«}   {F}dV+j{5u}   {<D}^  +  X{8«},r{p}, 

*  &  '"'  Eq.2.  1 

=  j({te}T{o}  +  {l>u}Tp{u}  +  {SuycAu}}iV 

Ve 

where  {8m}  and  {8s}  are  respectively  small  arbitrary  displacements  and  their 
corresponding  strains,  {F}  are  body  forces,  J®}  are  prescribed  surface  tractions  (which 
are  typically  nonzero  over  only  a  portion  of  surface  Se),  {/?},  are  concentrated  loads  that 
act  at  total  of  n  points  on  the  element,  {8m}  ,  is  the  displacement  of  the  point  at  which 
load  {/?},  is  applied,  p  is  the  mass  density  of  the  material,  Cd  is  a  material-damping 
parameter  analogous  to  viscosity,  and  the  volume  integration  is  carried  out  over  the 
element  volume  Ve. 

Using  usual  Finite  Element  notation,  we  may  write  the  continuous  displacement 
field  {«},  which  is  a  function  of  both  space  and  time,  and  its  first  two  time  derivatives, 
as 

\u}=[N]{d}  {&}-[Np}  {ii}=[N]{d}  Eq.2. 2 

In  Eqs.  2.2  the  so  called  shape  functions  [N]  are  functions  of  space  only,  and  the 
nodal  DOF  {d}  are  functions  of  time  only.  Thus  Eqs.  2.2  represent  a  local  separation  of 
variables.  Combination  of  Eqs.  2.1  and  2.2  yields 


MT 


\[B]T{cj}dV+  \p[NY[N]dv{d}+  \kd[N)T[N]dv{d) 

Eq.  2.  3 

0 


-j[N]'{F}dV-j[N]T{®}dS-Z{p 

Ve  Se  H 

in  which  it  has  been  assumed  that  the  concentrated  loads  {p},  are  applied  only  at  the 
nodal  points  locations.  Since  {6d\  is  arbitrary,  Eq.  2.3  can  be  written  as 
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[m]{d}  +  [c]{d}  +  {r'"<}={r<"}  Eq.  2.  4 

where  the  element  mass  and  damping  matrices  are  defined  as 

[m]=  \P[N]T[N]dV  Eq.2.5 

V 

\c]=\cd[N]T\N]dV  Eq.2.6 

V 

and  the  element  internal  forces  and  external  loads  vector  are  defined  as 

{r*}  =l[B]'{c}dV  Eq.2.7 

Ve 

{r"'}=l[N{{F}dV+  l[N]'{^}dS  +  fj{p};  Eq.2.8 

Ve  Se  '=1 

Equation  2.4  is  a  system  of  coupled,  second-order,  ordinary  differential 
equations  in  time  and  is  called  a  finite  element  semidiscretization  because  although  the 
displacements,  {a?},  are  discrete  functions  of  space,  they  are  still  continuous  functions  of 
time.  Methods  of  dynamic  analysis  focus  on  how  to  solve  this  equation.  Modal  methods 
focus  on  how  to  uncouple  the  equations,  transforming  the  NDOF  coupled  system  into  A' 
uncoupled  SDOF  systems,  each  of  which  can  be  solved  independently  of  others.  More 
details  about  this  method  of  analysis  can  be  found  in  Chapter  4.  Direct  integration 
methods  discretize  Eq.  2.4  in  time  to  obtain  a  sequence  of  algebraic  equations. 

Structure  matrices  [M],  [C],  and  {/?'"'}  are  constructed  by  standard  Finite 
Element  Method  procedures,  i.e.  conceptual  expansion  of  element  matrices  [m],  [c],  and 
{r""}  to  "structure  size"  followed  by  addition  of  overlapping  coefficients,  in  the  same 
way  it  is  done  for  assembling  stiffness  matrices  in  static  problems.  However,  the  exact 
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manner  in  which   {R""}   is  computed  depends  on  the  dynamic  analysis  procedure 
adopted. 

When  Eqs.  2.5  and  2.6  are  evaluated  using  the  same  shape  functions  [N]  as  used 
in  the  displacement  field  interpolation  (Eqs.  2.2),  the  results  are  called  consistent  mass 
and  consistent  damping  matrices.  These  matrices  are  symmetric.  On  the  element  level, 
they  are  generally  full,  but  on  the  structure  level,  they  have  the  same  sparse  form  as  the 
structure  stiffness  matrix.  When  p  and  Cd  are  nonzero,  consistent  matrices  [m]  and  [c] 
are   positive    definite.    Using    the   mass    matrix    for   example,    the   kinetic    energy 

l{</}  [m]{t/|  is  positive  definite  for  any  nonzero  {d}. 

In  typical  structural  analysis  we  are  more  interested  in  dry  fiction  and  hysteresis 
loss,  than  in  viscous  damping.  It  is  still  not  well  understood  how  the  damping 
mechanisms  develop  in  structures,  so  from  a  practical  standpoint  Eq.  2.6  does  not 
correctly  represent  structural  damping. 

The  internal  force  vector,  Eq.  2.7,  represents  loads  at  nodes  caused  by  straining 
of  material.  Equations  2.4  and  2.7  are  valid  for  both  linear  and  nonlinear  material 
behavior;  that  is,  in  Eq.  2.7,  {a}  could  be  a  nonlinear  function  of  strain  or  strain  rate. 
For  linearly  elastic  material  behavior,  {a}  =  [E][B]{d]  and  Eq.  2.7  becomes 

{'"'}=  MM  Eq.2.9 

where  the  usual  definition  of  the  stiffness  matrix  holds,  that  is, 

[*]=  \[B]T[E\B]dV  Eq.2.10 

Ve 

when  Eq.  2.10  is  used,  Eq.  2.4  becomes 
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[m]{d}  +  [c]{d}+[k]{d}  =  {r'"}  Eq.2.11 

which  can  be  interpreted  as  saying  that  external  loads  are  equilibrated  by  a  combination 
of  inertial,  damping,  and  elastic  forces.  For  the  assembled  structure,  from  Eq.  2.11  we 
get  the  equation  of  motion  for  linear  systems, 

[M]{D}  +  [C]{ />}+[£]{  D}  =  {R"'}  Eq.2.12 

where  {Rext}  corresponds  to  loads  {R}  of  a  static  problem,  but  is  in  general  a  function  of 
time.  Or,  returning  to  Eq.  2.4,  equations  of  the  assembled  structure  can  be  written  in  the 
alternative  form 

[M]{ D\  +  [C]{ D\  +  { RiM }  =  { R"' }  Eq.  2.  13 

which  does  not  require  that  the  material  be  linearly  elastic  and  represents  the  equation  of 
motion  for  nonlinear  systems. 

Equations  of  Motion  for  Ground  Motion 

It  is  now  opportune  to  derive  the  equations  of  motion  for  structural  systems 
subjected  to  ground  motion.  Consider  the  tower  shown  in  Fig.  2-1,  modeled  as  a 
cantilever  beam  with  concentrated  masses  at  the  nodes 

The  displacement  of  the  ground  is  denoted  by  Dg,  the  total  (or  absolute) 
displacement  of  the  mass  ms  by  Efs  and  the  relative  displacement  between  this  mass  and 
the  ground  by  Dj.  At  each  instant  of  time  these  displacements  are  related  by 

D'J(t)  =  Dj(l)  +  Dg(t)  Eq.2.14 

Such  equations  for  all  the  n  masses  can  be  combined  in  vector  form: 
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{D}'(t)  =  {D}(t)  +  {D}g(t)[l] 
where  [1]  is  a  vector  of  order  n  with  each  element  equal  to  unity. 


Rigid-body 
Motion 


Eq.  2. 15 


Fig.  2-1 .  Tower  subjected  to  ground  motion  -  after  Chopra  (1995) 


Only  the  relative  motion  [D]  between  the  masses  and  the  base  due  to  structural 
deformations  produce  elastic  and  damping  forces  (i.e.  the  rigidy  body  component  of  the 
displacement  of  the  structure  produces  no  internal  forces).  Thus  Eq.  2.13  is  still  valid, 
however  the  inertia  forces  are  related  to  the  total  acceleration  {/)'}  and  from  Eq.  2.15 
we  can  write 


and  substituting  this  value  back  into  Eq.  2.13  we  obtain 
[M]{D}+[C]{D}  +  {R-}  =  {PtJf} 


Eq.  2.  16 


Eq.  2.  17 
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The  external  force  vector  now  becomes  the  effective  earthquake  forces  vector, 
and  is  given  by 


Mr 


[MljoJ  Eq.2.18 


A  generalization  of  the  preceding  formulation  is  useful  if  all  the  DOF's  of  the 
system  are  not  in  the  direction  of  the  ground  motion,  or  if  the  earthquake  excitation  is 
not  identical  at  all  the  structural  supports  (see  Chapter  5  for  more  details).  In  this  general 
approach  the  total  displacement  of  each  mass  is  expressed  as  its  displacement  Lfj  due  to 
static  application  of  the  ground  motion  plus  the  dynamic  displacement  Dj  relative  to  the 
quasi-static  displacement: 

{D}'(t)  =  {D}  +  {D}'(t)  Eq.2.19 

The  quasi-static  displacements  can  be  expressed  as  {/)}*(/)  =  ({d}  (r),  where 
the  influence  vector  I  represents  the  displacements  of  the  masses  resulting  from  static 
application  of  a  unit  ground  displacement;  thus  Eq.  2.19  becomes 

{£>}'(/)  =  {D}(t)  +  e{D}f(t)  Eq.2.20 

The  equations  of  motion  are  obtained  as  before,  except  that  Eq.  2.20  is  used 
instead  of  Eq.  2.15,  resulting  in 

[Mp}+  [C]{D}+  {r™<  }=  -{M^{b\  Eq.  2.  21 

Now  the  effective  earthquake  forces  are 

WMwMAJ  Eq.2.22 
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To  help  illustrate  the  concept,  consider  the  inverted  L-shaped  frame  with  lumped 
masses  subjected  to  horizontal  ground  motion  shown  in  Fig.  2-2.  Assuming  the  elements 
to  be  axially  rigid,  the  three  DOF's  are  as  shown.  Static  application  of  Dg=  1  results  in 
the  displacements  shown  in  Fig.  2-2.  Thus  £={11  0}T  in  Eq.  2.21,  and  Eq.  2.22 
becomes 


P.H  (')  =  -\M\Df  (0  =  -Dg  (0  m2  +  m3 

0 


Eq.  2.  23 


"': 


'"i 


TO, 


o- 


(^ 


12=1 


£ 


I3=0 


Il=l 


E555S3 

W 


Mw 


-(mj+wj^rr) 


Stationary  base 


a) 


b) 


c) 


Fig.  2-2.  Support  motion  of  an  L-shaped  frame. 

a)  L-shaped  frame;  b)  influence  vector  £  :  static  displacements  due  to  Dg=l;  c) 
effective  load  vector  -  after  Chopra  (1995) 


Note  that  the  mass  corresponding  to  D2  =  1  is  toj+TOj,  because  both  masses  will 
undergo  the  same  acceleration  since  the  connecting  beam  is  axially  rigid.  The  effective 
forces  in  Eq.  2.23  are  shown  in  Fig  2-2.  Observe  that  the  effective  forces  are  zero  in  the 
vertical  DOF's  because  the  ground  motion  is  horizontal. 


18 
Mass  Matrices.  Consistent  and  Lumped 

A  mass  matrix  is  a  discrete  representation  of  a  continuous  distribution  of  mass. 
A    consistent    element    mass    matrix    is    defined    by    Eqs.    2.5    -    that    is,    by 

[m]  =  Jp[A^]  [A'JrfK.  It  is  termed  "consistent"  because  [N]  represents  the  same  shape 

v 

functions  as  are  used  in  the  displacement  field  interpolation,  and  to  generate  the  element 
stiffness  matrix.  A  simpler  formulation  is  the  lumped  mass  matrix,  which  is  obtained  by 

placing  particle  masses  m\  at  nodes  /  of  an  element,  such  that  J]  m,  is  the  total  element 
mass.  Particle  "lumps"  have  no  rotary  inertia  unless  rotary  inertia  is  arbitrarily  assigned, 
as  is  sometimes  done  for  the  rotational  DOF  of  beams  and  plates.  A  lumped  mass  matrix 
is  diagonal  but  a  consistent  mass  matrix  is  not.  The  two  formulations  have  different 
merits,  and  various  considerations  enter  into  deciding  which  one,  or  what  combination 
of  them,  is  best  suited  to  a  particular  analysis  procedure.  The  mass  matrix  for  a  3D 
uniform  beam  element,  which  is  used  to  model  the  pier  and  piles,  and  for  a  shell 
element,  which  is  used  to  model  the  pile's  cap,  are  developed  next. 
Mass  Matrix  for  the  Uniform  3D-Beam  Element 

Consistent 

The  formulation  found  here  is  given  by  Przemieniecki  (1968).  As  a  local 
coordinate  system,  consider  the  system  shown  in  Fig.2.3.  The  origin  is  at  node  1  with 
the  ox  axis  taken  along  the  length  of  the  beam  and  with  the  oy  and  oz  axis  as  the 
principal  axes  of  the  beam  cross  section.  The  matrix  N  for  this  element  consists  of 
twelve  displacements,  six  deflections  and  six  rotations,  that  is, 


19 


d={dl,d2,...,du 


Eq.  2.  24 
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Fig.  2-3.  3D  Beam  element 


Using  the  engineering  theory  of  bending  and  torsion  and  neglecting  shear 
deformations,  we  can  easily  show  that  the  matrix  N  in  the  relationship  {u}=[N]{d}  is 
given  by  Eqs.  2.25. 
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The  nondimensional  parameters  used  in  these  equations  are 


5=J 


Eq.  2.  26 


where  L  is  the  length  of  the  beam.  The  matrix  N  can  then  be  substituted  into  Eq.  2.5,  and 
performed  integration  over  the  whole  volume  of  the  element.  The  resulting  12  x  12 
consistent  mass  matrix  is  given  by 
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Eq.  2.  27 


where  the  matrix  terms  with  the  moments  of  inertia  Iy  or  I:  represent  rotatory  inertia  and 
the  terms  with  the  polar  moment  of  inertia  Jx  represent  torsional  inertia  of  the  element. 
Lumped 

The  lumped  mass  matrix  for  the  uniform  beam  segment  of  a  three-dimensional 
frame  element  is  simply  a  diagonal  matrix  in  which  the  coefficients  corresponding  to 
translatory  and  torsional  displacements  are  equal  to  one-half  of  the  total  inertia  of  the 
beam  segment  while  coefficients  corresponding  to  flexural  rotations  are  assumed  to  be 
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zero.  The  diagonal  lumped  mass  matrix  for  the  uniform  beam  of  distributed  mass 
m  =  pA  and  polar  mass  moment  Im  =  pJx  of  inertia  per  unit  of  length  may  be  written 
conveniently  as 


mi 


[l   1   1    Imlm   0   0   1   1   1   Iu/m   o   o] 


Eq.  2.  28 


Tables  2-1  and  2-2  compare  the  accuracy  of  the  consistent  and  lumped  mass 
formulations  using  finite  elements. 


Table  2-1.  Natural  frequencies  of  a  uniform  cantilever  beam:  Consistent-Mass  Finite 
Element  and  exact  solution 


Number  of  Finite  Elements,  Ne 

Mode 

1 

2 

3 

4 

5 

Exact 

1 

3.53273 

3.51772 

3.51637 

3.51613 

3.51606 

3.51602 

2 

34.8069 

22.2215 

22.1069 

22.0602 

22.0455 

22.0345 

3 

75.1571 

62.4659 

62.1749 

61.9188 

61.6972 

4 

218.138 

140.671 

122.657 

122.320 

120.902 

5 

264.743 

228.137 

203.020 

199.860 

Source:  Chopra  (1995). 

Table  2-2.  Natural  frequencies  of  a  uniform  cantilever  beam:  Lumped-Mass  Finite 
Element  and  Exact  Solution 


Number  of  Finite  Elements,  Ne 

Mode 

1 

2 

3 

4 

5 

Exact 

1 

2.44949 

3.15623 

3.34568 

3.41804 

3.45266 

3.51602 

2 

16.2580 

18.8859 

20.0904 

20.7335 

22.0345 

3 

47.0294 

53.2017 

55.9529 

61.6972 

4 

92.7302 

104.436 

120.902 

5 

153.017 

199.860 

Source:  Chopra  (1995). 
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Mass  Matrix  for  the  Shell  Element 

Because  in  FLPIER  the  pile's  cap  is  modeled  as  true  rectangular  9-node  shell 
elements,  we  will  limit  the  formulation  to  this  particular  type  of  element. 
Consistent 

Consider  the  true  rectangular  9-node  shell  element  shown  in  Fig.  2-4  below. 


Fig.  2-4.  True  9-node  rectangular  element 


Now  consider  the  following  mapping 
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Fig.  2-5.  Mapping  for  a  true  rectangular  9-node  shell  element 


It  is  easy  to  verify  that  the  shape  functions  N  for  each  node  are 
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tf,=-(|l-lXl|-l)HTl 


A,s=-^(H  +  1)(M-1)(ti-1)ti 


iV2=-(jl  +  lXTl-l)l«l 


N(,  =  --(v+m)(n-m+v 


JV,=-(h  +  1Xti +  Dl»1 


Af7=--(n  +  iXn-iXii  +  i>n  +  i 


^=-(^-1x11  +  1)1*11 


A's  =--(^)(h-1)(ti+1Xti-1) 


Af9=(M  +  l)(n-l)(ri+lXTl-l) 


Eqs.  2.  29 


Centerline 


Fig.  2-6.  Shell  element  of  uniform  thickness 


If  now  we  consider  for  each  node  the  six  DOF's  illustrated  in  Fig.  2-6,  which 
define  the  shell  element,  the  matrix  of  shape  functions  [N],  recalling  again  the 
relationship  {»}=  [N]{d],  takes  the  form 
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N,...N. 


N,...N, 


N-- 


N....N. 


N,...N. 


Eq.  2.  30 


N,...N0 


N,...N„ 


Note  that  [A/]'[A/]  has  dimensions  54  x  54,  which  are  the  correct  dimensions  for 
the  mass  matrix  considering  the  shell  element  illustrated  in  Fig.  2-6. 

It  is  easily  verified  that  for  the  true  rectangular  element  illustrated  in  Fig.  2-4  the 
change  of  coordinates  from  u  and  T)  to  x  and  y  can  be  expressed  as 


x  =  a\i  y  =  br\ 

and  recalling  for  convenience  Eq,  2.5 


[m]=  jp[N]T[N]dV 


Eq.  2.  31 


Eq.  2.  32 


If  we  now  consider  uniform  the  thickness  1  and  the  mass  density  p  over  the 
entire  element,  the  differential  volume  dVcan  be  written  dV=tdxdy  and  Eq.  2.5  now 
becomes 


[m]  =  pt\[N\'[N\dxdy 

v 

the  differencials  dx  and  dy  can  obtained  directly  from  Eq.  2.3 1, 

dx  =  ad\i  dy  =  bdr\ 

and  Eq.  2.32  can  be  rewritten 

i  i 
{m}  =  ptab\\[N]'[N}d\idr\ 


Eq.  2.  33 


Eq.  2.  34 


Eq.  2.  35 
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The  integral  in  Eq.  2.35  can  be  easily  evaluated  by  means  of  Gaussian 
Quadrature  (the  reader  is  referred  to  Appendix  B  for  more  details  on  this  procedure), 
therefore  any  element  of  the  mass  matrix  can  be  obtained  using  Eqs.  2.30  and  Eqs.  2.35 
and  this  completes  the  formulation  of  the  mass  matrix  for  the  true  rectangular  9-node 
shell  element. 
Lumped 

Lumping  the  mass  for  a  beam  element  is  a  process  that  seems  to  be  possible  by 
intuition  and  physical  insight,  however,  for  higher-order  elements,  like  the  shell 
element,  or  elements  of  irregular  shape,  intuition  can  be  risky.  Accordingly,  systematic 
schemes  for  lumping  are  necessary.  In  FLPIER  the  HRZ  scheme  is  used. 

The  HRZ  scheme  (Cook  et  al.,  1989)  is  an  effective  method  for  producing  a 
diagonal  mass  matrix.  It  can  be  recommended  for  arbitrary  elements.  The  idea  is  to  use 
only  the  diagonal  elements  of  the  consistent  mass  matrix,  but  to  scale  then  in  such  a  way 
that  the  total  mass  of  the  element  is  preserved.  Specifically,  the  procedural  steps  are  as 
follows  (Cook  et  al.,  1989): 

1)  Compute  only  the  diagonal  coefficients  of  the  consistent  mass  matrix. 

2)  Compute  the  total  mass  of  the  element,  m. 

3)  Compute  a  number  s  by  adding  the  diagonal  coefficients  ma  associated  with 

the  translational  DOF  (but  not  rotational  DOF,  if  any)  that  are  mutually 
parallel  and  in  the  same  direction. 

4)  Scale  all  the  diagonal  coefficients  by  multiplying  them  by  the  ratio  m/s,  thus 

preserving  the  total  mass  of  the  element. 
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Following  this  procedure  for  the  true  rectangular  9-node  shell  element  with 
uniform  thickness  /,  results  in  the  following  lumped  mass  distribution 
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Fig.  2-7.  Lumped  mass  matrix  at  the  nodes  of  true  rectangular  9-node  shell  element. 
Numbers  shown  are  fractions  of  the  total  element  mass  at  each  node 


Remarks  about  the  mass  matrix 

Cook  et  al.  (1989)  makes  important  remarks  about  the  mass  matrix.  The  first  one 
is  the  fact  that  the  mass  matrix  chosen  must  correctly  represent  the  mass  distribution  on 
the  element,  because  the  product  [m]{rf}  or  [M]{/3}  must  yield  the  correct  total  force  on 
an  element  according  to  Newton's  law  F  =  ma  when  {d}  represents  a  rigid-body 
translational  acceleration.  The  second  remark  is  about  the  consistent  and  lumped  mass 
matrix.  While  the  consistent  mass  matrix  is  always  positive  definite,  the  same  can  not  be 
said  about  the  lumped  mass  matrix.  If  it  contains  zeros  or  negative  entries  in  the 
diagonal,  then  it  is  called  positive  semi-definite.  This  may  or  may  not  cause  some  matrix 
operations  to  give  strange  results.  He  also  suggests  that  a  consistent  mass  matrix  may  be 
more  suitable  for  flexural  problems,  while  lumped  mass  matrices  usually  yield  natural 
frequencies  that  are  less  than  the  exact  values. 
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The  main  reason  for  the  development  of  lumped  mass  schemes  was  to  save 
memory  space  in  computational  calculations.  In  the  past  this  was  a  problem,  but  today 
computers  are  much  faster  and  memory  has  become  cheap  and  abundant.  Therefore  the 
use  of  consistent  mass  formulation  is  justified,  since  it  avoids  certain  instabilities  in  the 
matrices  operations. 

Damping 

Damping  in  structures  is  not  viscous,  i.e.  is  not  proportional  to  velocity;  rather  it 
is  due  to  mechanisms  such  as  hysteresis  in  the  material  and  slip  in  connections.  These 
mechanisms  are  not  yet  well  understood.  Moreover,  these  mechanisms  are  either  too 
difficult  to  incorporate  into  the  analysis,  or  they  make  the  equations  computationally  too 
expensive.  Therefore  with  the  actual  limited  knowledge  about  damping  mechanisms, 
viscous  damping  is  usually  adopted  in  most  analysis.  Comparisons  of  theory  and 
experiment  show  that  this  approach  is  sufficiently  accurate  in  most  cases. 

Damping  in  structures  can  be  considered  in  two  ways: 

1)  phenomenological  damping  methods,  in  which  the  actual  physical  dissipative 

mechanisms  such  as  elastic-plastic  hysteresis  loss,  structural  joint  friction,  or 
material  microcracking  are  modeled. 

2)  spectra!  damping  methods,  in  which  viscous  damping  is  introduced  by  means 

of  specified  fractions  of  critical  damping  (Critical  damping,  for  which  the 

damping   ratio   is   t,  =  1 ,  marks   the   transition   between   oscillatory   and 

nonoscillatory  response). 

The   first   class   of  methods   requires   detailed    models    for   the   dissipative 

mechanisms  and  almost  always  result  in  nonlinear  analyses.  In  the  second  class  of 
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methods,  experimental  observations  of  the  vibratory  response  of  structures  are  used  to 
assign  a  fraction  of  critical  damping  as  a  function  of  frequency,  or  more  commonly,  a 
single  damping  fraction  for  the  entire  frequency  range  of  a  structure.  The  damping  ratio 
\  depends  on  the  material  properties  at  the  stress  level.  For  example,  in  steel  piping 
%  ranges  from  about  0.5%  at  low  stress  levels  to  about  5%  at  high  stress  levels.  In  bolted 
or  riveted  steel  structures,  and  in  reinforced  or  prestressed  concrete,  \  has  the 
approximate  range  2%  to  15%  (Cook  et  al.,  1989). 

A  popular  spectral  damping  scheme,  called  Rayleigh  or  proportional  damping  is 
used  to  form  the  damping  matrix  [C]  as  a  linear  combination  of  the  stiffness  and  mass 
matrices  of  the  system,  that  is 

[c]  =  a[A:]+p[M]  Eq.2.36 

where  a  and  (5  are  called,  respectively,  the  stiffness  and  mass  proportional  damping 
constants.  Matrix  [C]  given  by  Eq.  2.36  is  an  orthogonal  damping  matrix  because  it 
permits  modes  to  be  uncoupled  by  eigenvectors  associated  with  the  undamped 
eingenproblem.  The  relationship  between  a,  p,  and  the  fraction  of  critical  damping  ^  is 
given  by  the  following  equation. 

If P 


§  =  2la0,+»J  E"-2-37 


damping  constants  a  and  p  are  determined  by  choosing  the  fractions  of  critical  damping 
(%\  and  tp)  at  two  different  frequencies  (to,  and  co2)  and  solving  simultaneous  equations 
for  a  and  p.  Thus: 

a=2(^<B2-%co,)/((o22-(0l2)  Eq.2.38 
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P  =  2co,co,(^co2  -^co,)  /(co22  -co,2)  Eq.  2.  39 

The  damping  factor  a  applied  to  the  stiffness  matrix  [K]  increases  with 
increasing  frequency,  whereas  the  damping  factor  p  applied  to  the  mass  matrix  [M] 
increases  with  decreasing  frequency.  For  structures  that  may  have  rigid-body  motion,  it 
is  important  that  the  mass-proportional  damping  not  be  excessive. 

Usually,  the  natural  frequencies  coi  and  k>2  are  chosen  to  bound  the  design 
spectrum.  Therefore  coi  is  taken  as  the  lowest  natural  frequency  of  the  structure,  and  C02 
is  the  maximum  frequency  of  interest  in  the  loading  or  response.  Cook  et  al.  (1989) 
suggests  a  value  of  30  Hz  as  the  upper  frequency  for  seismic  analysis,  because  the 
spectral  content  of  seismic  design  spectra  are  insignificant  above  that  frequency. 

Estimating  Modal  Damping  Ratios 

Because  damping  is  still  an  unknown  subjected  the  estimate  of  damping  rations 
still  presents  some  challenge.  Recommended  damping  values  are  given  in  Table  2-3  for 
two  levels  of  motion:  working  stress  levels  or  stress  levels  no  more  than  one-half  the 
yield  point,  and  stresses  at  or  just  below  the  yield  point.  For  each  stress  level,  a  range  of 
damping  values  is  given;  the  higher  values  of  damping  are  to  be  used  for  ordinary 
structures,  and  the  lower  values  are  for  special  structures  to  be  designed  more 
conservatively.  In  addition  to  Table  2-3,  recommended  damping  values  are  3%  for 
unreinforced  masonry  structures  and  7%  for  reinforced  masonry  construction.  It  is 
important  to  note  that  the  recommended  damping  ratios  given  in  Table  2-3  can  only  be 
used  for  the  linearly  elastic  analysis  of  structures  with  classical  damping.  This  implies 
that  all  the  material  components  of  the  structure  are  still  behaving  in  their  linear-elastic 
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phase.  For  structures  subjected  to  strong  motions,  that  will  lead  to  crushing  of  concrete 
or  yielding  of  steel,  characterizing  nonlinear  material  behavior,  hysteretic  damping  must 
be  added  to  the  analysis  through  nonlinear  force-deformation  relationships. 

Table  2-3.  Recommended  damping  rations  for  structures 


Stress  Level 

Type  and  Condition 
of  structure 

Damping  Ratio(%) 

Working    stress,    no    more 
than  about  Vi  yield  point 

Welded    steel,    prestressed 
concrete,       well-reinforced 
concrete        (only        slight 
cracking) 

2-3 

Reinforced    concrete    with 
considerable  craking 

3-5 

Bolted  and/or  riveted  steel, 
wood  structures  with  nailed 
or  bolted  joints 

5-7 

At  or  just  below  yield  point 

Welded    steel,    prestressed 
concrete  (without  complete 
loss  in  prestress) 

5-7 

Prestressed    concrete    with 
no  prestress  left 

7-10 

Reinforced  concrete 

7-10 

Bolted  and/or  riveted  steel, 
wood  structures  with  bolted 
joints 

10-15 

Wood  structures  with  nalied 
joints 

15-20 

Source:  Chopra  (1995) 


Mass  Condensation 


A  useful  tool  to  decrease  the  number  of  DOF's  in  the  static  analysis  of  a  system 
without  losing  accuracy  is  static  condensation.  In  this  approach  some  degrees  of 
freedom  of  the  structure  are  chosen  as  master  degrees  of  freedom  and  the  remaining 
ones  are  called  slave  degrees  of  freedom.  The  choice  for  the  master  DOF  is  concerned 
with  those  DOF  that  give  a  better  representation  of  the  system,  i.e.  in  a  building  where 
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the  beams  are  much  stiffer  than  the  columns  the  shear  DOF's  would  be  a  good  choice, 

instead  of  the  rotations.  The  condensed  stiffness  matrix  obtained  includes  the  effects  of 

the  slave  degrees  of  freedom,  which  can  be  recovered  at  any  time  during  analysis.  The 

same  approach  can  be  applied  to  the  mass  matrix,  however  dynamic  condensation 

(Miller,  (1981),  Paz  (1985))  is  not  exact,  as  will  be  shown  later  in  this  section,  but  it  can 

give  good  results  if  some  rules  are  observed  in  the  modeling.  The  approach  described 

here  is  given  by  Meirovitch  (1980,  pp.371-372): 

Let  us  write  the  equations  for  the  potential  and  kinetic  energy  for  a  system  in  the 
matrix  form: 


V =^[D]T[K][D] 


T  =  -[D]T[M][D] 


Eq.  2.  40 
Eq.  2.  41 


and  divide  the  displacement  vector  [D]  into  the  master  displacement  vector  q2 
and  slave  displacement  vector  qi,  or 


[D]- 


Eq.  2.  42 


Then   the   stiffness   matrix    [K]    and   mass   matrix    [M]    can   be   partitioned 
accordingly,  with  the  result 


Eq.  2.  43 
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K„ 
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M„ 

w12p, 

Mn 

M22{q2 

Eq.  2.  44 


where  K21-K12   and  M2i=Mi2T.  The  condition  that  there  be  no  applied  forces  in 
the  direction  of  the  slave  displacements  can  be  written  symbolically  in  the  form 


dV 
dq, 


-q,Ku+q2K2l 


Eq.  2.  45 
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where  the  equation  implies  equilibrium  in  the  direction  of  the  slave 
displacements.  Solving  Eq.  2.45  for  qi,  we  obtain 

9i  =-KuKnl2 


Eq.  2.  46 


which  can  be  used  to  eliminate  qi  from  the  problem  formulation. 

Equation  2.46  can  be  regarded  as  a  constraint  equation,  so  that  the  complete 

displacement  vector  q  can  be  expressed  in  terms  of  the  master  vector  q2  in  the 

form 


q~Pqt 


Eq.  2.  47 


where  P  is  a  rectangular  constraint  matrix  having  the  form 

/ 


Eq.  2.  48 


in  which  /  is  a  unit  matrix  of  the  same  order  as  the  dimension  of  q2.  Introducing 
Eq.  2.47  into  Eqs.  2.40  and  2.41,  we  obtain 


r  =  2«iK2«2 


T  =  \qlM2q2 


Eq.  2.  49 
Eq.  2.  50 


where  the  reduced  stiffness  and  mass  matrices  are  simply 


Kt  —  P  KP  -  K22  -  KnKu  K2] 


Eq.  2.  51 


M,  =PTMP=  M22  -K2\K^'M2I  -  MnKfiKx+XaXnKiXnKzi  Eq-2.52 


The  matrix  Mi  is  generally  known  as  the  condensed  mass  matrix. 

What  is  being  sacrificed  as  a  result  of  the  condensation  process? 

To  answer  this  question,  let  us  consider  the  complete  eigenvalue  problem,  which 

can  be  separated  into 


Kn9,  +Knq2  =k(Muq,  +  Ml2q2) 
Kix9i+K22q2  =\(M2lq,  +  M22q2) 


Eq.  2.  53 
Eq.  2.  54 


Solving  Eq.  2.54  for  qi,  we  have 

q,  =(XM2l-K2l)(K22-XM22)q2 


Eq.  2.  55 
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so  that,  introducing  Eq.  2.55  into  Eq.  2.53,  we  obtain 

(jC„  -K,2K2lK2l)q2  =  x(M22  -KnKu'M2l  -  Ml2KjK2,  +  KnKn' MuKu'K2l)q,  + 

Eq.  2.  56 

Examining  Eqs.  2.56  and  2.52  we  conclude  that  the  condensation  used  earlier  is 
ignoring  second  and  higher  order  terms  in  X  in  Eq.  2.56,  which  can  be  justified  if 
the  coefficients  of  X2,  Xs,  ...  are  significantly  smaller  than  the  coefficient  of  X. 
For  this  to  be  true  we  must  have  the  entries  of  Mn  and  Mn  much  smaller  than 
the  entries  of  Kn  and  K}j.  Physically,  this  implies  that  the  slave  displacements 
should  be  chosen  from  areas  of  high  stiffness  and  low  mass.  Moreover  the  nodes 
that  carry  a  time-varying  load  should  be  retained  as  master. 

FLIPIER  condenses  the  stiffness  and  mass  to  the  top  of  the  piles.  While  this 
procedure  is  exact  for  the  stiffness,  note  that  it  is  not  for  the  mass.  This  is  because  in 
pier  structures  the  slave  DOF's  are  located  in  areas  of  high  mass  concentration,  like  the 
pile's  cap.  Also  note  that  the  mass  of  the  superstructure  on  the  pier,  is  usually  modeled 
as  lumped  masses  at  the  top  of  the  pier.  In  the  case  of  earthquake  loading  for  example 
note  that  we  have  the  loading  function  acting  on  slave  DOF's,  what  is  not  acceptable  in 
this  approach.  Because  of  these  limitations  for  mass  condensation,  this  approach  is  not 
recommended  for  the  dynamic  analysis  of  bridge  piers.  FLPIER  was  then  modified  to 
allow  a  "full"  analysis  of  the  structure,  where  neither  the  stiffness  nor  the  mass  matrices 
are  condensed  to  the  top  of  the  piles.  A  typical  full  and  the  respective  condensed  version 
of  a  structure  are  shown  in  Fig.  2-8. 

In  a  typical  condensed  static  analysis,  the  stiffness  and  loads  of  the 
superstructure  are  condensed  to  the  top  of  the  piles.  The  condensed  analysis  is  then 
carried  out.  At  the  end  of  the  analysis  the  superstructure's  forces  and  displacements  are 
recovered  and  the  analysis  is  terminated. 


34 


Condensed  masses 


Cap 


\ 


Condensed 
structure 


Fig.  2.8.  Full  and  condensed  versions  of  the  structure 


Time-History  Analysis.  Direct  Integration  Methods 


In  direct  integration  methods  or  step-by-step  methods,  a  finite  difference 
approximation  is  used  to  replace  the  time  derivatives  appearing  in  Eq.  2.12  or  2.13  at 
various  instants  of  time.  Direct  integration  is  an  alternative  to  modal  analysis  methods. 
For  many  structural  dynamics  and  wave  propagation  problems,  including  those  with 
complicated  nonlinearities,  direct  integration  is  easier  to  implement.  In  direct 
integration,  the  approach  is  to  write  the  equation  of  motion  (2.12),  at  a  specific  instant  in 
time, 


[M]{D}n+[C]{D}ii+[K]{D}„  =  {R'"}n 


Eq.  2.  57 


where  the  subscript  n  denotes  time  nAt  and  At  is  the  size  of  the  time  increment  or  time 
step.  The  absence  of  time  step  subscripts  on  matrices  [M],  [C],  and  [K]  in  Eq.  2.57 
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implies  linearity.  For  problems  with  material  or  geometric  nonlinearity,  [K]  is  a  function 
of  displacement  and  therefore  of  time  as  well.  Accordingly,  from  Eq.  2.57, 

[M]{zJ}n  +[C]{£)} <  +  {jf*^  =  {R'"}i:  Eq.  2.  58 

{/?""}„  is  the  internal  force  vector  at  time  n  At  due  to  straining  of  material.  It  is  obtained 
by  assembling  element  internal  force  vectors,  {/""}„,  given  by  Eq.  2.7  using  {o}„.  For 
nonlinear  problems,  {R'"'}„  is  a  nonlinear  function  of  {£)}„  and  possibly  time  derivatives 
of  {D}„,  if  the  strain  rate  is  an  issue.  For  linear  problems,  the  internal  force  vector  is 
given  by  the  relationship  {R'"'}„=[K]{D}„.  In  Eq.  2.58  [M]  and  [C]  are  taken  as  time- 
independent,  although  for  some  problems  these  may  also  be  nonlinear.  In  this  work  [M\ 
and  [C]  remain  constant  during  the  analysis,  and  the  internal  force  vector,  {/?""},  is  only 
a  function  of  the  displacements  {D}„.  The  nonlinear  relationship  for  {a}„  will  be 
developed  later. 

Different  methods  for  direct  integration  of  Eqs.  2.57  and  2.58  can  be  categorized 
as  explicit  or  implicit.  The  first  category,  the  explicit  methods,  has  the  form 

{0}„+1  =f({D}n,{D}jD}jD}  „_„...)  Eq.2.59 

and  hence  permit  {D}„+i  to  be  determined  in  terms  of  completely  historical  information 
consisting  of  displacements  and  time  derivatives  of  displacements  at  time  n  At  and 
before.  The  main  characteristic  of  explicit  methods  is  the  fact  that  the  next 
approximation  for  the  displacements  is  based  only  on  the  known  previous 
approximations  for  the  displacements  and  their  respective  derivatives.  Note  the  use  of 
the  equilibrium  condition  at  time  n.  The  Central  Difference  is  an  example  of  an  explicit 
method. 
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The  second  category,  the  implicit  methods,  has  the  form 

{«},„,  =/({i>}„tl,{o}/i+|,{D}„,...)  Eq.2.60 

and  hence  computation  of  {D}„+i  requires  knowledge  of  the  time  derivatives  of  {D}„+i, 
which  are  unknown.  The  main  characteristic  of  the  implicit  methods  is  the  fact  that  the 
next  approximation  for  the  displacements  depends  on  unknown  values  of  their 
derivatives.  Note  that  the  equilibrium  condition  is  used  at  time  step  n+1.  Newmark's 
and  the  Wilson-Theta  are  examples  of  implicit  methods. 

There  is  vast  literature  about  the  advantages  and  disadvantages  of  each  approach, 
the  reader  is  referred  to  Cook  et  al.  (1989),  Chopra  (1995),  Paz  (1985),  Craig  (1981),  or 
Crandall  (1956),  for  a  more  extensive  discussion  on  these  approaches.  Generally 
speaking  under  certain  conditions  the  implicit  methods  are  more  stable  than  the  explicit 
methods.  Because  an  implicit  method  was  used  in  this  research  we  will  limit  our 
discussion  to  this  class  of  methods. 

Numerical  Evaluaton  of  Dynamic  Response.  Newmark's  Method 

As  mentioned  earlier  implicit  methods  are  those  in  which  the  approximation  for 
the  next  displacements  {D}„+i  depends  on  unknown  values  for  its  time  derivatives.  The 
main  advantage  of  the  implicit  methods  is  the  fact  that  most  of  the  useful  methods  are 
unconditionally  stable  and  have  no  restriction  on  the  time  step  size  other  than  required 
for  accuracy. 

In  1959  N.  M.  Newmark  developed  a  family  of  time-stepping  methods  based 
assumptions  for  the  variation  of  the  acceleration  over  the  time  step.  The  first  method  is 
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called  average  acceleration  and  is  shown  in  Fig.  2-9.  Successive  application  of  the 
Trapezoidal  Rule  leads  to  Eqs.  2.61  to  2.65. 


it   ■■ 

"m 

"i 

| — ► 

/ 

At 
Fig.  2-9.  Average  acceleration 


«CO«-(fiU,+fi,) 


w(t)  =  «,  +  -  («„,+",) 


A/ 


m(t)  =  u{  +  u^t  +  ~r(w/+i  +  «,) 


At' 
uM  =  u:  +  ulAt  +  ——(u^l  +U;) 


Eq.  2.  61 
Eq.  2.  62 
Eq.  2.  63 
Eq.  2.  64 
Eq.  2.  65 


The  second  method  is  called  linear  acceleration  and  is  illustrated  in  Fig.2-10. 
The  Trapezoidal  Rule  is  also  used  to  obtain  Eqs.  2.66  to  2.70. 


Km 


At 
Fig.  2-10.  Linear  acceleration 
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"(t)  =  ",+—(»,>,  +",) 


m(t)  =  ii,  +m,t  H (it..  +ii|) 


A/ 

",♦1  =»,+y  ("/.i+"/) 


X2         T3 

m(t)  =  uf  +  in  +  ii,  — -  + — (uM  -  u, ) 


«.+|  =  «.  +u;A/  +  (A<)     —"/+1+-«( 
Vo  3 


Eq.  2.  66 
Eq.  2.  67 
Eq.  2.  68 
Eq.  2.  69 
Eq.  2.  70 


The  Newmark's  family  of  methods  can  be  summarized  into  the  following  two 
equations: 

",+,  ="i+[0-y)"/+y",>iW 


Uui  =  u.  +U.A/  + 


|-P 


A/; 


Eq.  2.  71 
Eqs.  2.  72 


where  y  and  p  are  parameters  that  can  be  determined  to  obtain  integration  accuracy  and 
stability.  Note  that  when  y  =  14  and  P  =  1/6,  this  method  reduces  to  the  linear 
acceleration  method.  Newmark  originally  proposed  as  an  unconditionally  stable  scheme, 
the  constant  average  acceleration  method,  in  which  case  y  =  'A  and  p  =  %.  This  can  be 
shown  considering  that  Newmark's  method  is  stable  if  (Chopra,  1995): 


At         1  1 

•  <- 


T„     7i  V2  Vy-2P 
For  y  =  Vi  and  p  =  %  this  condition  becomes 


Eq.  2.  73 
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A/ 

—  <  oo  Eq.  2.  74 

T. 


This  implies  that  the  average  acceleration  method  is  stable  for  any  Af,  no  matter 
how  large,  as  mentioned  before;  however,  it  is  accurate  only  if  Af  is  small  enough.  For  y 
=  Vt  and  p  =1/6,  Eq.  2.73  indicates  that  the  linear  acceleration  method  is  stable  if 


—  <  0.551  Eq.2.  75 

T. 


For  the  solution  of  displacements,  velocities  and  acceleration  at  time  i+\  we 
now  consider  the  equilibrium  equation  also  at  time  i+l: 

mM  +  CuM  +  KuM  =  FM  Eq.  2.  76 

Solving  Eq.  2.72  for  «/+1in  terms  of  w,+i  and  substituting  in  Eq.  2.71,  equations 
for  uMand  uM  in  terms  of  the  unknown  displacements  «,+/  only  are  obtained. 
Substituting  these  equations  into  Eq.  2.76,  a  system  of  equations  is  obtained,  which  can 
be  solved  to  obtain  u,+a,: 

(b„M  +  btC  +  K)u,.  =  FM  +  M(b0ut  +  b2U,  +  b3ui)+  Cfau,  +  bAu,  +  65H,)Eq.  2.  77 

where 


— Urlb,  *—;b,  =—;b,=- — \;b,  =^--\;b.  =A/M  Eq.  2.  78 

pA/2     '      PA/     2      PA/     '      2p  4      p  5 


and  finally  all  the  quantities  at  time  (+1  can  be  written  as 

"i+i  =  fy>(",>i  ~"/)  +  *7"/  +*8",  Eq.  2.  79 
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"m  =  "/  +  b,iit  +  bloiiltl  Eq.  2.  80 

uM  m  u,  +  («w  -  u, )  Eq.  2.  81 


where 


b6  =  b0 ;  A7  =  -62 ;  *8  =  -b} ;  69  =  A/(7-y;  ;  bw  =  Y- A'  E1-  2-  82 

Equation  2.77  can  be  now  written  in  the  condensed  form: 

KuM=p(t)  Eq.2.83 

where  the  effective  stiffness  K  is  given  by 

K  =  b0M  +  btC  +  K  Eq.2.84 

and  the  effective  load  vector  p(t)  is 

p(t)  =  FM  +  M(baut  +b2u,  +  63u, )+  C(6,«;  +*4»(  +  Mi )  E<1-  2-  85 

which  completes  the  formulation  for  linear  systems. 

Choice  of  time  step  At 

The  unconditional  stability  of  the  average  acceleration  method  may  lead  some 
analysts  to  adopt  larger  time  steps  because  of  economic  needs,  implying  a  solution  that 
may  not  be  accurate,  because  unconditional  stability  does  not  mean  unconditional 
accuracy.  Cook  et  at.  (1989)  suggests  the  following  expression  for  the  time  step  At: 

A/ <(2h/g>„)/20  =  0.3/co,,  Eq.  2.  86 

where  co„  is  the  highest  frequency  of  interest  in  the  loading  or  response  of  the  structure. 
However,  he  suggests  that  in  the  case  of  convergence  difficulties,  the  analysis  should  be 
repeated  with  a  smaller  time  step  for  additional  assurance  of  a  correct  solution. 
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Nonlinear  Problems 

In  structural  analysis,  a  problem  is  nonlinear  if  the  stiffness  matrix  or  the  load 
vector  depends  on  the  displacements.  The  nonlinearity  in  structures  can  be  classified  as 
material  nonlinearity  (associated  with  changes  in  material  properties,  as  in  plasticity)  or 
as  geometric  nonlinearity  (associated  with  changes  in  configuration,  as  in  large 
deflections  of  a  slender  elastic  beam).  In  general,  for  a  time-independent  problem 
symbolized  as  [K]{D}={R],  in  linear  analysis  both  [K]  and  {R}  are  regarded  as 
independent  of  {£>},  whereas  in  nonlinear  analysis  {K}  and/or  {R}  are  regarded  as 
functions  of  {D}.  In  this  dissertation  we  will  address  the  nonlinearities  associated  with 
changes  in  material  properties. 
Analysis  of  the  Nonlinear  Response  using  Newmark's  Method 

In  this  section  the  notation  for  SDOF  systems  is  used  to  simplify  the  approach, 
the  extension  to  MDOF  systems  is  done  later,  although  the  concept  is  exactly  the  same. 
In  the  numerical  evaluation  of  the  response  of  a  dynamic  system  we  go  from  time  step  ;', 
where  the  equation  of  motion  can  be  written 

mw,  +  cu,  +  (fs)t  =  />,  Eq.  2.  87 

to  time  step  ;'  + 1 ,  where  the  dynamic  equilibrium  can  be  written 

m"„.  +  c",+i  +  0)w  =  pM  Eq- 2-  88 

where  (fs),  is  the  system's  resisting  force  at  time  i.  For  a  linear  system  (fs),  =  hi,. 
However  for  a  nonlinear  system  the  resisiting  force  fs,  would  depend  on  the  prior  history 
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of  displacements  and  velocities.  It  is  then  necessary  to  consider  the  incremental 
equilibrium  equation,  the  difference  between  Eqs.  2.84  and  2.85,  can  be  written: 

mAiif  +  cAu,  +  (Afs)t  =  Ap,  Eq.  2.  89 

The  incremental  resisting  force  can  be  written 

(A/S), ,=(*,),«,  An,  Eq.2.90 

where  the  secant  stiffness  (£,)sec,  shown  in  Fig.  2-11,  cannot  be  determined  because  «,+i 
is  not  known.  If  however  we  make  the  assumption  that  over  a  small  time  step  At  the 
secant  stiffness  (£,)sec,  can  be  replaced  by  the  tangent  stiffness  (£,)tan,  then  Eq.  2.90  can 
be  rewritten: 

(&fs)i={kl)TAul  Eq.2.91 

The  incremental  dynamic  equilibrium  equation  is  now: 

mAii,  +  cAu,  +  (k,)T  Am,  =  Ap,  Eq.  2.  92 

Equation  2.92  suggests  that  the  analysis  of  nonlinear  systems  can  be  done  by 
simply  replacing  the  stiffness  matrix  k  by  the  tangent  stiffness  (ki)Tlo  be  evaluated  at  the 
beginning  of  each  time  step.  However  this  procedure  for  constant  time  steps  At  can  lead 
to  unacceptable  results  for  two  reasons: 

a)  The  tangent  stiffness  was  used  instead  of  the  secant  stiffness. 

b)  The  use  of  a  constant  time  step  delays  detection  of  the  transitions  in  the  force- 

deformation  relationship. 
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Fig.  2-11.  Secant  and  Tangent  approaches.  After  Chopra(1995) 

These  errors  can  be  minimized  by  using  an  iterative  procedure  within  each  time 
step.  The  idea  is  to  guarantee  dynamic  equilibrium  before  going  to  the  next  time  step. 
We  must  then  solve  the  equation 


**,♦,  =  PC) 


where  now  the  effective  stiffness  k,  becomes 


Eq.  2.  93 


k.  =  k:  +  b„m  +  b,c 


Eq.  2.  94 


For  convenience  we  drop  the  subscript  i  in  k;  and  replace  it  by  T  to  emphasize 
that  this  is  the  tangent  stiffness,  and  we  can  rewrite  Eq.  2.94  as 


kT  =  kT  +  b0m  +  blc 


Eq.  2.  95 


The  effective  force  p(t)  is  now  given  by 

P(')  =  /('),+,  +  m(b1ui  +  bjU, )  +  c(b,u,  +  b5u,tl )  -  r, 
where  r,  is  the  internal  force  at  time  step  i. 


Eq.  2.  96 
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Since  we  are  using  the  tangent  estimate  for  the  stiffness,  we  must  iterate  within 
each  time  step  to  guarantee  dynamic  equilibrium.  The  first  step  is  to  apply  the  effective 
force  p(t)  and  get  the  first  approximation  for  the  next  displacement  u  (+;.  Associated 
with  this  displacement  is  the  new  tangent  stiffness  Icr  and  the  true  force  f,  which 
includes  the  internal  force  r,  the  inertia  force//,  and  the  damping  force  fc  at  iteration 
k.  An  out-of  balance  force  defined  as  Ap(t)  +1  =  p(t)  -/is  then  generated.  The  new 
effective    stiffness    for    the    system    is    computed    and    the    system    of   equations 

kk  Au*  =  Ap*  is  solved  to  obtain  the  incremental  displacements  A«* ,  which  is  added  to 
the  previous  displacement  u,,  giving  a  new  estimate  for  the  displacements,  velocities  and 
accelerations  at  time  /+1.  The  dynamic  equilibrium  is  checked  again  with  the  new 
values  and  another  out-of  balance  force  is  generated.  The  process  continues  until  the 
out-of-balance  force  is  within  a  specified  tolerance.  Figure  2-12  helps  to  understand  the 
process.  This  proccess  is  known  as  the  Newton-Raphson  method. 

If  the  reader  is  familiar  with  numerical  methods,  will  notice  that  what  is  defined 
as  the  The  Newton-Raphson  method  in  some  references  (Cook  et  al.  (1989),  Chopra 
(1995)  and  Bathe  (1996))  and  above  is  in  fact  a  Newton-Raphson  approach.  Although 
irrelevant  to  this  discussion  the  author  found  important  to  define  the  Newton-Raphson 
method  from  a  mathematical  point  of  view. 

The  Newton-Raphson  method  is  a  mathematical  method  for  finding  roots.  The 
method  relies  on  the  fact  that  the  function  fix),  and  its  first  derivative  f(x),  and  second 
derivative  /'  '(x)  are  continuous  near  a  root  p.  Note  how  rigorous  the  mathematical 
definition  for  the  method  is,  since  this  was  not  assumed  in  the  previous  application. 
With  this  information  it  is  possible  to  develop  algorithms  that  produce  sequences  {**} 
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that  converge  faster  to  p.  Note  however  \ha\.f(x)  and  its  first  and  second  derivatives  must 
be  continuous  for  the  method  to  work.  The  convergence  rate  of  the  Newton-Raphson 
method  is  quadratic.  This  method  is  also  called  the  tangent  method  because  the  slope  of 
the  curve  is  used  in  the  formulation.  The  recursive  formula  for  the  method  is: 


*■(**!  )  =  **-4rr4. 


fix*) 


Eq.  2.  97 


and  the  sequence  will  converge  top. 


Fig.  2-12.  Newton-Raphson  Method 


However  sometimes  it  may  be  difficult  if  not  impossible  to  implement  Newton's 
method  if  the  first  derivative  is  complicated.  The  secant  method  may  then  be  used.  The 
secant  method  is  the  same  as  the  Newton's  method  except  that  the  first  derivative  is 
replaced  by  its  approximation:  the  slope  of  the  line  through  the  two  previous  points. 
Note  that  when  successive  points  get  close,  the  method  becomes  unstable  because 
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division  by  zero  may  occur.  It  can  be  shown  that  the  convergence  rate  is  approximately 
1.62,  being  slower  than  Newton's  method.  The  recursive  formula  for  the  secant  method 
is  shown  in  Eq.  2.97.  Note  the  clear  connection  between  the  mathematical  and 
engineering  terms  used  to  define  both  approaches.  More  detailed  information  about  both 
methods  can  be  found  in  Mathews  (1987). 

nx^x,-l&L  (fa)  /ex*.-*-)     Eq.2.98 

Xn  ~Xn-l 

The  extension  to  MDOF  systems  is  immediate  by  replacing  all  the  scalar 
quantities  by  its  respective  vector  equivalents.  Note  however  that  for  a  SDOF  system  the 
tangent  stiffness  kr  is  easily  obtained.  Such  a  simple  evaluation  is  not  available  if  there 
are  MDOF.  However,  in  practice,  as  it  will  be  discussed  later,  the  physics  of  the 
problem  allows  us  to  calculate  the  tangent-stiffness  matrix  [K,].  In  a  MDOF  context  N-R 
iteration  involves  repeated  solution  of  the  equations  [X,],{AD},+/={Ai?},+;,  where  the 
tangent  stiffness  matrix  [K,]  and  load  imbalance  {AS}  are  updated  after  each  cycle.  The 
solution  process  seeks  to  reduce  the  load  imbalance,  and  consequently  {AD},  to  zero. 
The  internal  force  in  the  equation  of  motion,  Eq.  2.35,  is  written  as 

k"'L+,  ={*'"'  I  +  MaD}  Eq.2.99 

where 

{AD}  =  {D}„tl  -  JD}„  Eq.  2.  1 00 

Combining  Eqs.  2.99  and  2.100  with  the  equations  of  motion,  Eq.  2.35,  and  the 
trapezoidal  rule,  we  obtain 


47 


where 


and 


[^]{AD}={^}„+,  Eq.2.101 


[K'ff\  =  A.[M]+L[C]+[Ki]  Eq.2.102 


Note  that  [£",]  must  be  predicted  using  {£>}„  (and  possibly  {£)}  if  strain  rates 
effects  are  important)  and  must  be  factored  at  least  once  each  time  step  during  nonlinear 
response.  If  [Ay  is  not  an  accurate  prediction  of  the  true  tangent-stiffness  matrix  from 
time  nAt  to  time  (n  +  l)Ar,  then  the  solution  of  Eq.  2.101  for  {AD}  will  be  in  error.  The 
error  in  nodal  forces,  that  is  the  residual,  is  given  by  the  imbalance  in  the  equation  of 
motion  as 

{ir}.  {*-},„,  -[m](4+i  -[cpU  -{/r'L,  Eq-  2.  mm 

where  the  internal  force  vector  {R'"'}„+i  is  computed  element-by  element.  The  process 
stops  when  the  out-of-balance  force  vector  {S"T}  is  smaller  than  a  specified  tolerance. 

Finally,  considering  all  the  topics  discussed  in  this  section,  the  procedure  for 
nonlinear  dynamic  analysis  can  be  summarized  in  the  algorithm  described  next. 
Nonlinear  Dynamic  Analysis  Algorithm 

Step  (1).  Form  the  initial  stiffness  matrix  K,  mass  matrix  M,  and  damping  C  for 
the  structure. 

Step  (2).  Compute  the  constants  for  the  numerical  method  chosen: 
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'ilson-9 

Newmark 

h               6        .3,            3 

°~(QAtf  "~9A, 
62  =  2&i;Z>3  =  2 

6. 

62 

..     1     .b  =_§_ 
cxAr2  '  '     a  At 
=l/oAf,63  =  l/2a-l 

64  =2;6,  =  9Af/2 
bt=bJQ;b7  =-62/G 

65  =l-3/9;6,  =  Ar/2 
610=A/2/6 

ft. 
6, 

ftin 

=  --l;65=A/(--2)/2 
a                    a 

=  60;A,  =  -62 

=  -ft3;ft,  =  Af(l-8) 

=  8A/ 

For  the  Newmark' s  Method  -  Average  Acceleration:  a  =  V*  and  5  =  !4. 

For  the  Newmark's  Method  -  Linear  Acceleration:  a  =  1/6  and  5  =  'A. 

For  the  Wilson-9  Method:  9  =  1 .4  (9  =  1 .0  for  Newmark's). 

Step  (3).  Initialize  ii0,u0  and  u0- 

Step  (4).  Form  the  effective  stiffness  matrix  using  the  initial  stiffness  matrix  K: 

K  =  bt)M  +  blC  +  K 

Step  (5).  Beginning  of  time  step  loop. 

Step  (6).  Form  the  effective  load  vector  for  the  current  time  step: 
£«*  =  F,  +6 (Fltil  -F,)  +  M(b2u,  +  b,u, )  +  Cib.ii,  +  b5u, ) - R, 
Step  (7).  Solve  for  the  displacement  increments 

5u  =  [k]'fm 

Step  (8).  Beginning  of  dynamic  equilibrium  loop. 

Step  (9).  Evaluate  approximations  for  ii,u  and  u: 

K3m  =b08u'~'  -b2ii,  -b,ii, 


"nQAI 


=  bt5u'  '  -b.it,  -b5u, 


"',~J>&,  =  "„8u'"' 
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Step  (10).  Evaluate  actual  tangent  stiffness  K1'1  and  internal  forces 

R'~'  for  all  the  elements  in  the  structure. 

Step  (11).  Evaluate  new  effective  stiffness  matrix  K,  based  on 

actual  tangent  stiffness  K!~' . 

Step  (12).  Evaluate  the  out-of-balance  dynamic  forces: 

X~L  =  f,  +o(fm  -f,)-[mu'-l  +o&L  + jcU 

Step  (13).  Evaluate  the  ith  corrected  displacement  increment: 

a«' =  [£;-'}"  V,^ 

Step  (14).  Evaluate  the  corrected  displacement  increments: 
8m'  =  &u'~'  +  Au' 

Step  (15).  Check  for  the  convergence  of  the  iteration  process: 
NO 


Return  to  Step  (8) 


YES 


Step  (16).  Return  to  Step  (5)  to  process  the  next  time  step. 

In  this  dissertation  the  full  Newton-Raphson  Method  is  applied.  In  this  approach 
the  tangent  stiffness  matrix   kr   is  updated  for  every  iteration,  in  contrast  with  the 

Modified  Newton-Raphson  Method  in  which  the  tangent  stiffness  kT  is  update  once 
every  time  step  (it  remains  constant  during  the  iterations  performed  within  each  time 
step).  This  improves  convergence,  but  additional  computational  effort  is  required  in 
forming  a  new  tangent  stiffness  matrix  kT  and  factorizing  it  at  each  iteration  cycle.  The 
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iterations  within  each  time  step  proceed  until  the  value  of  the  out-of-balance  force 
Sf^j,  is  smaller  than  a  specified  tolerance  t.  The  mass  matrix  M  and  damping  matrix  C 
are  considered  to  be  constant  throughout  the  analysis. 

As  mentioned  earlier  in  this  Chapter  the  correct  choice  of  time  step,  which 
depends  on  the  lowest  natural  frequency  of  interest  for  the  system,  enforces  stability  and 
accuracy  to  the  analysis.  Note  however  that  for  a  nonlinear  analysis,  the  case  is  that  due 
to  yielding  of  steel  or  crushing  and  cracking  of  concrete,  a  lower  stiffness  will  be 
obtained  during  the  analysis,  which  will  result  in  even  lower  values  for  the  lowest 
natural  frequency  of  interest.  The  reader  should  have  this  in  mind  when  choosing  an 
adequate  time  step  for  a  nonlinear  dynamic  analysis. 


CHAPTER  3 
DISCRETE  ELEMENT  MODEL  AND  MATERIAL  HYSTERESIS 


Discrete  Element  Derivation 

A  representation  of  the  discrete-element  used  in  FLPIER  is  shown  in  Fig.  3-1.  It 
was  developed  by  Mitchell  (1973)  and  modified  by  Andrade  (1994).  The  center  bar  can 
both  twist  and  extent  but  is  otherwise  rigid.  The  center  bar  is  connected  by  two  universal 
joints  to  two  rigid  and  blocks.  The  universal  joints  permit  bending  at  the  quarter  points 
about  the  y  and  z  axes.  Discrete  deformational  angle  changes  ¥/,  4^,  ^3  and  Tj  occur 
corresponding  to  the  bending  moments  Mi,  M2,  M3  and  M4,  respectively.  A  discrete 
axial  shortening  (8)  corresponds  to  the  axial  thrust  (T),  and  the  torsional  angle  ¥5 
corresponds  to  the  torsional  moments  in  the  center  bar  M5. 
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Fig.  3-1 .  Representation  of  discrete  element.  After  Hoit  et  al.,  1 996 
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52 
Element  Deformation  Relations 

In  Fig.  3-2,  W1-W3  and  w7-w9  represent  the  displacements  in  the  x-,  y-  and  z- 
directions  at  the  left  and  right  ends,  respectively.  The  displacements  w4  and  ww 
represent  axial  twists  (twists  about  the  x-axis)  at  the  left  and  right  ends,  respectively;  and 
finally  w5-w6  and  wn-wl2  represent  the  angles  at  the  left  and  right  end  blocks  about  the  x 
and  z  axes,  respectively.  Based  on  a  small  displacement  theory  we  can  write: 

h 
«  =  w3-w9--(w5  +  wn)  Eq.3. 1 

h 
s=ws-w2--(w6  +  wl2)  Eq.3. 2 

The  elongation  of  the  center  section  of  the  element  is  calculated  as  follows: 

8  =  wi  ~  wi  Eq.  3.  3 

The  angle  changes  for  the  center  section  about  the  z  and  y  axes  are  the  defined  as 

Eq.  3.  4 

Eq.  3.  5 

The  discretized  vertical  and  horizontal  angle  changes  at  the  two  universal  joints 
are  then 

Ti=e.-H'6  Eq.3. 6 

^2  =  w3  -62  Eq.  3.  7 

li'3  =  wi2-ei  Eq.3. 8 


s 

ws-w2 

w6  +  wa 

h 

h 

2 

s 

h 

2 
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Eq.  3.  9 

and  the  twist  in  the  center  part  of  the  element  is  defined  as 

^5  =  w]0  -  w4  Eq.  3. 10 

Therefore,  the  internal  deformations  of  the  discrete  element  model  are  uniquely 
defined  for  any  combination  of  element  end  displacements.  The  curvature  for  small 
displacements  at  the  left  and  right  universal  joints  about  the  y  and  z  axes  are  defined  as 
follows: 

At  the  left  joint 

<t»i  ='*',//»  Eq.3. 11 

h='¥2/h  Eq.3. 12 

At  the  right  joint 

^=%/h  Eq.3. 13 

<|>4=4'4//i  Eq.3. 14 

The  axial  strain  at  the  center  of  the  section  is  given  by 

z<=Th  E<»-315 
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Fig.  3-2.  Discrete  element  displacements.  After  Hoit  et  al.  1996 


Integration  of  Stresses  for  Nonlinear  Materials 


For  a  beam  subjected  to  both  bending  and  axial  loads,  it  is  assumed  that  the 
strains  vary  linearly  over  the  area  of  the  cross  section.  This  assumption  enables  the 
strain  components  due  to  bending  about  the  z  and  .y-axes,  and  the  axial  strain  to  be 
combined  using  super  position.  Examples  of  these  three  components  are  represented 
separately  in  Fig.  3-3(a-c)  and  combined  in  Fig.  3-3(d)  also  shown  in  Fig.  3-3(d)  is  a 
differential  force,  dFt,  acting  on  a  differential  area,  dAj,  relationship  for  the  material 

dF,=  0,(14,  Eq.3. 16 

Finally,  Fig.  3-3(e)  represents  the  stress-strain  curve.  Then,  for  the  left  joint,  the 
relationship  for  the  strain  at  any  point  in  the  cross  section  is 
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Then  to  satisfy  equilibrium 
M,  =  J  \dFt.Y,  =  J  jofidA, 

A  A 

My  =  JJrff;.^  =  \\o,Z,dA, 

.4  /I 


Eq.3.  17 

Eq.  3.  18 
Eq.  3. 19 
Eq.  3.  20 


a)  Strain  due  to 
z-axis  bending 


e)  Stress-strain  relationship 


-^ 


b)  Strain  due  to 
/■axis  bending 


^p^ 


t^^r 


c)  Strain  due  to 

axial  thrust 


d)   Combined  strains 


Fig.  3-3.  Various  components  of  total  strain  in  the  section.  After  Hoit  et  al.,  1996 


Numerical  integration  of  Eqs.  3.18,  3.19  and  3.20  is  done  using  Gaussian 
Quadrature.  To  use  the  method  of  Gaussian  Quadrature,  the  function  being  integrated 
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must  be  evaluated  at  those  points  specified  by  the  position  factor.  These  factors  are  then 
multiplied  by  the  appropriate  weighting  factors  and  the  products  accumulated.  Fig.  3-4 
shows  a  square  section  with  25  integration  points  (a  5  x  5  mesh).  The  actual  number  of 
defaults  integration  points  for  a  square  section  is  set  at  81  (a  9  x  9  mesh).  For  a  steel  H- 
section  the  default  number  of  points  is  60.  For  circular  sections,  the  section  is  divided 
into  circular  sectors  (12  radial  divisions  and  five  circumferencial  divisions  as  shown  in 
Fig.  3-5),  totaling  60  points.  The  sections  are  integrated  at  the  centroids  of  each  sector 
using  weighting  factors  of  1 .0.  The  stress  in  all  steel  bars  is  evaluated  at  the  centroid  and 
a  weighting  factor  of  1  is  used  for  each  bar.  When  a  circular  void  is  encountered  in  a 
square  section,  the  force  is  first  computed  on  the  unvoided  section  and  then  the  force 
that  would  be  acting  on  the  voided  circular  area  is  computed  and  subtracted  from  the 
force  computed  for  the  unvoided  section.  Circular  sections  with  voids  are  divided  into 
sectors  omitting  the  voided  portion.  This  method  of  dividing  the  sections  into  points, 
getting  the  strains  and  stresses  at  each  point  and  then  integrating  to  get  forces  is  usually 
called  fiber  modeling. 
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Fig.  3-4.  Rectangular  section  with  integration  points.  After  Hoit  et  al.,  1996 
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Steel  bar  (lxl  Integration) 


Concrete  subdivided 
into  sectors 


Note:  Integration  po  nts  { lxl )  for  concrete  are 
at  die  geometric  centroids  of  each  sector 


Fig.  3-5.  Circular  section  with  integration  points.  After  Hoit  et  al,  1996 


Even  for  a  nonlinear  material  analysis,  the  torsional  moment  M  is  assumed  to 
be  a  linear  function  of  the  angle  of  twist  S?  ,  and  the  torsional  stiffness  GJ ,  where  J  is 
the  cross  section  torsional  constant;  and  G  is  the  material  shear  modulus,  resulting  in 
the  following  expression  for  the  torsional  moment 

M  =  GJVI1h  Eq.3.21 

In  this  discrete  approach  the  curvature  is  evaluated  at  two  deformational  joints 
inside  the  element.  These  points  are  located  at  the  quarter  points  from  each  node  along 
the  element  length.  Therefore  the  effective  position  of  any  plastic  hinge  that  might  form 
in  the  structure  is  restricted  to  these  two  locations.  This  should  not  cause  any  practical 
limitation  for  most  problems.  However,  it  should  be  considered  if  trying  to  match  a 
theoretical  solution  with  pure  plastic  hinges  at  theoretical  ends  of  members.  These 
discrete  joints,  at  which  the  deformations  are  concentrated,  correspond  to  the  integration 
points  along  the  length  of  the  element,  if  a  conventional  finite-element  solution  was 
being  made. 
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Element  End  Forces 

The  element  internal  forces  are  necessary  to  assemble  the  global  internal  force 
vector  necessary  for  equilibrium  of  the  equation  of  motion.  From  equilibrium  of  the 
center  bar 

r,-(j/4-J/,)/A  Eq.3.22 

V2=(M,-M,)/h  Eq.3.23 
and  from  equilibrium  of  the  end  bars 

f\  =  ~T  Eq.  3.  24 

fi  =  vi  Eq.  3.  25 

fi  =  ~vi  Eq.  3. 26 

f*=-M5  Eq.3.27 

/5=Ml+F2-A/2+r-A/2-wJ  Eq.3.28 

f6=M2+Vrh/2  +  T-h/2-w6  Eq.3.29 

/?  =  r  Eq.  3.  30 

/«  =  "Fi  Eq.  3.  31 

/»  =  F2  Eq.  3.  32 

/w  =  M5  Eq.  3.  33 

/„  =-M,  +  V2-h/2  +  T-h/2-wu  Eq.3.34 

fl2=M4+Vrh/2  +  T-h/2-wu  Eq.3.35 

where/;  -/5  and/7  -/9  are  the  acting  end  forces;  and/,  -/«  and/,,,  -fa  are  the  acting  end 
moments. 
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Element  Stiffness 

As  mentioned  in  Chapter  2  the  element  stiffness  may  change  for  either  nonlinear 
static  or  dynamic  analysis.  Therefore  its  is  necessary  to  evaluate  the  element  stiffness 
for  each  iteration  in  each  time-step.  The  procedure  adopted  uses  the  standard  definition 
of  the  stiffness  matrix;  for  an  element  having  n  DOF  the  stiffness  matrix  is  a  square 
matrix  [K]  of  dimensions  n  x  n,  in  which  Kti  is  the  force  necessary  in  the  rth  DOF  to 
produce  a  unit  deflection  of  they'th  DOF.  The  stiffness  computed  is  that  obtained  by  one 
of  the  two  methods  described  below.  The  transformation  of  the  discrete  element 
stiffness  matrix  to  global  coordinates  and  the  assembly  of  the  different  components  of 
the  global  stiffness  matrix  follow  standard  direct  stiffness  procedures. 

Secant  and  Tangent  Stiffness  of  the  Discrete  Element 

During  the  iteration  process,  the  element  stiffness  matrix  is  reevaluated  in  each 
new  deformed  position.  For  each  iteration,  the  stiffness  for  each  integration  point  along 
the  cross  section  within  an  element  is  stored.  Then,  on  12  subsequent  passes,  a  unit 
displacement  is  applied  to  each  element  DOF  -  keeping  all  other  DOF  fixed  -  and  the 
forces  corresponding  to  that  unit  displacement  are  calculated  over  the  cross  section  of 
the  element  as  described  earlier.  If  the  stiffness  for  each  of  the  integration  points  is 
defined  by  dividing  the  present  stress  by  the  present  strain  as  shown  in  Fig.  3-6,  then  it 
is  called  secant  stiffness.  On  the  other  hand,  if  the  stiffness  for  each  segment  is  defined 
by  the  slope  of  the  stress-strain  curve  at  the  specific  point  being  integrated,  then  it  is 
called  the  tangent  stiffness  as  also  illustrated  in  Fig.3-6.  Note  that  if  the  tangent  stiffness 
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is  used,  the  solution  of  the  nonlinear  problems  becomes  Newton-Raphson,  while  if  the 
secant  stiffness  is  used  we  have  a  secant  method  solution. 
Stress 


Strain 
Fig.  3-6.  Secant  and  tangent  material  stiffness 


Hysteresis  Models 

We  have  shown  so  far  how  it  is  possible  to  get  the  element  internal  forces  and 
form  the  updated  stiffness  for  each  iteration  in  the  previous  paragraphs.  But  we  have 
considered  only  the  case  in  which  the  element  is  subjected  to  loading.  For  dynamic 
analysis  usually  the  case  is  that  the  element  will  be  subjected  to  two  additional  phases: 
unloading  and  reloading.  It  is  important  to  notice  that  even  for  nonlinear  static  analysis 
these  two  new  phases  could  be  present  due  to  a  nonlinear  redistribution  offerees  on  the 
structure  that  could  cause  some  elements  to  be  over-loaded  and  others  to  be  under- 
loaded. The  curve  that  is  used  to  describe  this  behavior  is  called  a  hysteresis.  For  the 
inelastic  analysis,  a  proper  selection  of  hysteretic  models  for  the  materials  is  one  of  the 
critical  factors  in  successfully  predicting  the  dynamic  response  under  strong  motion. 
Several  models  have  been  proposed  in  the  past  for  reproducing  various  aspects  of 
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reinforced  concrete  behavior  under  inelastic  loading  reversals.  In  order  to  closely 
reproduce  the  hysteretic  behavior  of  various  components,  a  highly  versatile  model  is 
required  in  which  several  significant  aspects  of  hysteretic  loops  can  be  included,  i.e., 
stiffness  degradation,  strength  deterioration,  pinching  behavior  and  the  variability  of 
hysteresis  loop  areas  at  different  deformation  levels  under  repeated  loading  reversals. 
However,  the  model  should  also  be  as  simple  as  possible  since  a  large  number  of 
inelastic  spring  are  necessary  in  modeling  the  entire  structure,  and  additional  parameters 
to  describe  a  complicated  hysteresis  loop  shape  may  sometimes  require  excessive 
amount  of  information. 
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Fig.  3-7.  Models  for  hysteresis  loops  proposed  by  some  authors.  After  Mo,  1994 

Some  of  the  existing  popular  models:  Clough  (1966),  Fukada(1969),  Ayoama 

(1971),   Kustu   (1975),   Tani   (1973),   Takeda   (1970),   Park   (1984),   Iwan   (1973), 

Takayanagi  (1977),  Muto  (1973),  Atalay  (1975),  Nakata  (1978),  Blakeley  (1973),  and 
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Mo  (1988)  are  shown  in  Fig.3-7.  It  appears  that  most  of  the  available  models  are  aimed 
at  a  particular  type  of  component,  such  as  for  use  of  beams,  columns  or  shear  walls  only, 
and  therefore,  fall  short  of  the  versatility  required  for  modeling  practical  buildings 
having  a  large  number  of  different  components.  Most  of  these  models  were  obtained  by 
performing  curve  fits  to  experimental  data.  This  approach  results  in  really  good 
approximations  for  the  behavior  of  the  specific  member  being  studied,  but  lacks 
versatility  when  applied  to  a  real  structure. 

The  advantage  of  the  discrete  element  with  fiber  modeling  is  that  the  elements 
general  behavior  will  be  governed  by  its  material  properties,  instead  of  experimental 
observations,  making  this  procedure  useful  for  studying  multiple  cross  section 
configurations  under  more  general  load  histories. 

Material  Models 

Rather  than  trying  to  develop  a  new  element  that  would  model  this  behavior 
using  elasticity  and  plasticity  theories,  it  was  thought  to  change  the  discrete  element 
used  in  FLPIER  to  develop  such  behavior.  Based  on  various  studies  and  experiments 
(Chen  (1982),  Park  and  Paulay  (1975),  Mo  (1994),  Roufaiel  and  Meyer  (1987),  Park  et 
al.  (1972),  Agrawal  et  al.  (1965),  Sinha  et  al.  (1964-a),  Sinha  et  al.  (1964-b),  Kwak 
(1997),  Magdy  and  Mayer  (1985),  Soroushian  et  al.  (1986),  Nilsson  (1979),  Ozcebe 
(1989),  Peinzien  (1960),  and  Tseng  (1976))  about  the  uniaxial  cyclic  behavior  of 
concrete  and  steel,  and  the  behavior  of  R/C  members,  the  stress-strain  curves  of  these 
materials  were  modified  to  work  for  dynamic  analysis.  Shear  deformations  were  not 
included  because  it  was  not  included  in  the  works  cited.  Moreover  it  was  also  not 
included  in  the  formulation  of  the  discrete  element.  By  using  this  approach  it  was 
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possible  to  overcome  some  of  the  difficulties  found  in  the  various  hysteresis  models. 
Although  computationally  more  expensive,  this  fiber  modeling  approach  gave  the 
modified  discrete  element  the  versatility  not  found  in  many  models.  A  description  of 
these  modified  curves  follows  next. 

Uniaxial  Mild  Steel  Model 

Because  the  cyclic  behavior  of  steel  is  very  dependent  on  the  heat  from  which  it 
was  produced,  it  was  decided  that  rather  than  trying  to  predict  the  steel  behavior  based 
on  exponential  curves  (Agrawal  et  al.  (1965),  Shen  and  Dong  (1997)),  to  use  the  bilinear 
representation  for  the  steel  behavior.  This  represents  a  safe  lower  bound  solution  that  is 
adequate  for  most  of  the  construction  steel.  So  the  mild  steel  reinforcement  is  assumed 
to  be  perfectly  elastic-plastic  (no  hardening)  and  similar  in  both  tension  and 
compression  as  shown  in  Fig.  3-8.  The  parameters  needed  for  the  mild  steel  bilinear 
model  are  the  modulus  of  elasticity  Es  and  yield  stress /y.  The  rules  for  this  model  are  as 
follows  (refer  to  Fig.  3-8): 

Loading  is  represented  by  segment  a-b,  the  tangent  stiffness  is  the  initial 
modulus  of  elasticity  for  steel  Es  and  the  stress  is  given  by  a  =  Es.£  . 

Yielding:  If  e  >  Ey  yielding  occurs  and  is  represented  by  segment  b-c.  The 
tangent  stiffness  Es  is  equal  to  0  and  the  stress  a  =  f  ,  the  yielding  stress  for  steel.  The 
residual  strain  Er  is  given  by  £r  =s  -e   . 

Unloading  is  represented  by  segment  c-d,  the  tangent  stiffness  is  Es  and  the 
stress  is  given  by  <j  =  (e  -£,).£, . 
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Reloading  is  represented  by  segment  e-f,  the  tangent  stiffness  is  E,  and  the 
stress  is  given  by  o"  =  (e  -  e  r  ).£, . 

For  any  of  these  phases  the  secant  stiffness  Ess  is  given  by 


E=- 


Eq.  3.  36 


Fig.  3-8.  Elastic-perfectly  plastic  model  for  mild  steel 


Uniaxial  Monotonic  Concrete  Model  Used  in  FLPIER 

The  concrete  model  used  in  FLPIER  is  generated  based  on  the  values  of  the 
concrete  strength  fc  and  modulus  of  elasticity  of  concrete  Ec,  input  by  the  user.  The 
compression  portion  of  the  curve,  which  is  is  based  on  the  work  of  Wang  and  Reese 
(1993),  is  highly  non-linear  and  has  a  maximum  compressive  stress^ ",  which  is  related 
to  but  not  always  equal  to  the  compressive  strength  of  a  standard  test  cylinder,/,'.  Based 
on  experimental  research,  fc  is  taken  to  be  85%  of  fc\  the  maximum  cylinder 
compression  stress. 
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The  tension  side  of  the  curve  is  based  on  the  tension  stiffening  model  proposed 
by  Mitchell  (1973).  This  procedure  assumes  an  average  tensile  stress-strain  curve  for 
concrete.  The  stress  strain  relation  of  concrete  in  tension  is  very  close  to  linear  with 
cracking  occurring  at  a  small  rupture  stress/,,  The  high  stresses  actually  experienced  at 
tensile  cracks  in  the  concrete  will  not  be  reproduced  by  the  model.  However  the  average 
response  over  a  finite  length  of  beam  will  be  adequately  represented. 

Based  on  the  user  input  the  program  will  generate  the  concrete  curve  as  a  series 
of  points  connected  by  straight  lines  as  shown  in  Fig.  3-9.  Values  in  between  the  points 
are  obtained  by  interpolation  of  the  extreme  points  in  the  interval.  For  values  of  fc  =  6 
ksi  and  Ec  =  4615  ksi  the  strain  and  stress  values  for  the  concrete  stress-strain  curve 
generated  by  the  program  can  be  seen  in  Figs.  3-10  and  3-1 1  respectively. 


I13.in  25/  J/1999 


Fig.  3-9.  FLPIER  concrete  points 
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Fig.  3-10.  Concrete  strains 
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Fig.  3-11.  Concrete  stresses 


67 
Proposed  Models  for  the  Uniaxial  Inelastic  Cyclic  Behavior  of  Concrete 

Two  concrete  models  were  implemented  in  FLPIER.  One,  called  the  rational 
model,  is    based  on  the  work  of  Sinha  et  al.  (1964-a),  the  other  is  a  bilinear  model, 
similar  to  the  mild  steel  model  presented  earlier.  These  two  models  are  discussed  in 
more  details  next. 
Rational  Model 

Figure  3-9  shows  the  default  envelop  of  the  stress  strain  curve  supplied  by  the 
program,  which  is  a  function  of  /'cand  Ec  input  by  the  user.  This  curve  is  the  backbone 
of  the  model,  because  it  limits  the  value  of  the  stress  in  concrete  during  all  phases  in  the 
analysis.  The  compression  portion  of  the  concrete  curve  is  highly  nonlinear  and  is 
defined  by  the  Hognstead  parabola  up  to  a  stress  equal  to  0.85/'c.  Beyond  this  point,  a 
straight  line  is  adopted,  connecting  this  maximum  to  another  point  with  residual  stress 
of  0.20/V  and  strain  equal  to  4  £„,  as  shown  in  the  Fig.  3-12.  For  strains  greater  than  4 
£„,  the  stress  is  kept  at  a  constant  0.20/V  (residual  stress)  as  suggested  by  Chen  (1982). 
For  the  tension  portion  the  curve  is  assumed  linear  up  to  a  stress  of  /,  and  then  has  a 
tension-softening  portion  as  shown  in  Fig.  3-12.  The  tension-softening  portion  attempts 
to  account  for  the  uncracked  sections  between  cracks  where  the  concrete  still  carries 
some  stress.  The  value  of  /,  is  based  on  the  fixed  value  of  e,  shown  in  the  figure  and 
the  modulus  of  elasticity  Ecl  input  by  the  user.  For  English  units  this  will  give  a  value 
of  fr  =  15jfc  .  The  rules  for  the  rational  model  are  described  next. 
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Fig.  3-12.  Envelop  curve  for  concrete 

Loading 

Compression 

Compression    follows    the    curve    described    above,    a    typical    loading    in 
compression  phase  is  illustrated  in  Fig.  3-13.  The  equations  that  define  this  phase  are: 

When  £<e„  concrete  follows  the  Hongnestead  parabola  defined  by 


'-'4rMr 


Eq.  3.  37 


The  tangent  modulus  of  elasticity  is  defined  as: 
If e<E»  then 


Ec=Ect 
If  £>E„then 


Eq.  3.  38 
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£j=M| !_.!.]  Eq>3_39 


This  correction  was  necessary  because  the  derivative  of  the  Hongnestead 
parabola  gives  very  high  values  of  the  tangent  modulus  for  low  values  of  the 
compressive  strain  s,  which  caused  instabilities  to  the  model.  Based  on  research  (Chen 
(1982),  Park  and  Paulay  (1975))  it  was  found  that  at  about  30%  of/'c  it  is  reasonable  to 
assume  that  concrete  still  have  the  initial  tangent  modulus,  Eci.  Equation  3.39  is  just  the 
derivative  of  Eq.  3.37  with  respect  to  e. 

When  the  strain  e  >  e0,  the  concrete  enters  a  phase  called  softening.  In  this 
research  this  phase  is  defined  as: 

If  e0<e  <4  E„then 

„  (  Q-8/;' 


'     (4s,  -sj 

and  if  E>4  eu 


Eq.  3.  41 


/  =  0.20/c"  Eq.  3.  42 

and 

Ec  =  0  Eq.  3.  43 
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Fig.  3-13.  Typical  compression  loading 

Tension 

A  typical  loading  in  tension  is  shown  in  Fig.  3-14.  Tension  is  also  represented  by 
the  envelope  curve  described  earlier.  The  following  equations  define  loading  in  tension 

IfE<Er 


f  =  Ee 

Ec  =  Ecl 

Ife,  <£<£.,,/ then 


Eq.  3.  44 
Eq.  3.  45 


/  =  /,+ 


0-5/, 


(*-*,) 


Eq.  3.  46 
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E  = 


0.5/r 


Eq.  3.  47 


IfEs,/<8  <e„,  then 


/-Ir^-K-s.) 


Eq.  3.  48 


E  = 


as/. 


Eq.  3.  49 


and  finally  if  £  >  e„i  then 

/  =  0 
and 


Eq.  3.  50 


E  =0 


Eq.  3.  51 


stress 
(ksi) 


-0.0001     0    0.0001       0.0003       0.0005      0.0007       o.OO 

strain  in/in 


Fig.  3-14.  Typical  loading  in  tension 
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Unloading 

Typical  unloading  phases  in  tension  and  compression  are  illustrated  if  Fig.  3-15 
and  Fig.  3-16,  respectively.  The  expressions  for  compression  and  tension  unloading  are 
based  on  the  work  of  Sinha  et  al.  (1964-a)  defined  for  different  concrete  mixes.  Due  to 
the  lack  of  more  experimental  data,  it  is  assumed  that  for  stronger  concrete  (f'c  >  4  ksi), 
the  response  will  be  that  of  4  ksi  concrete.  If  new  experiments  are  done  for  stronger 
concrete  these  changes  can  be  easily  incorporated  into  the  program.  The  family  of 
unloading  curves  is  represented  by  second  order  curves.  From  experiments  for  various 
concrete  strengths  a  good  fit  was  obtained  with  the  formula: 

a  +  H  =  — (e  -  Xf  Eq.  3. 52 

where  H  and  J  are  experimental  constants  whose  values  for  the  three  mixes  used  are 
given  in  Table  3-1.  For  stress  in  units  of  ksi  and  strains  in  units  of  in/in,  and  X  is  a 
particular  parameter,  different  values  of  which  represent  different  members  of  the 
family.  To  determine  the  value  of  the  parameter  X,  for  a  curve  passing  though  a  certain 
point  in  the  stress-strain  plane  a,  ,E, ,  the  coordinate  values  are  substituted  into  Eq.  3.52 
which  is  then  solved  for  the  required  value  of  X,  leading  to  the  expression: 


X  =  E'+^-J(£'+^J-E'2 


The  tangent  modulus  Ec  for  unloading  is  taken  as  the  slope  for  two  consecutive 
points  in  the  unloading  curve: 


f.-f 

ijl—L-  Eq.  3.  54 
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Table  3-1.  Curves  coefficients 


f  c  (psi) 

H 

J 

K 

L 

3000 

0.07 

0.95 

3.42 

1.26 

3750 

0.09 

0.52 

2.52 

1.03 

4000 

0.10 

0.61 

4.61 

1.01 

There  are  two  choices  for  when  the  unloading  in  compression  crosses  the  strain 
axis.  The  first  option  assumes  that  a  gap  is  formed  in  compression,  and  concrete  will 
totally  unload  until  reloading  in  tension,  as  illustrated  in  Fig.  3-18.  The  second  option 
assumes  that  no  gap  is  formed,  and  concrete  will  go  straight  into  the  tension  reloading 
phase  from  compression,  as  shown  in  Fig.  3-17.  When  unloading  from  tension  a  gap  is 
always  formed,  and  concrete  will  not  go  into  compression  until  the  gap  is  closed,  as 
shown  in  Fig.  3-15. 
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Fig.  3-15.  Typical  unloading  in  tension 
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Fig.  3-16.  Typical  unloading  in  compression 
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Fig.  3-17.  Compression  unloading  with  gap 
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Fig.  3-18.  Compression  unloading  with  no  gap 


Reloading 

A  typical  reloading  in  compression  is  illustrated  in  Fig.  3-19.  The  reloading 
curves  in  either  compression  or  tension  are  represented  by  a  family  of  converging 
straight  lines;  accordingly,  an  expression: 

ct  +  K  =  7(6  +  L)  Eq.  3. 55 

was  chosen,  in  which  K  and  L  are  experimental  constants,  and  7,  is  a  parameter.  The 
value  of  7  can  be  found,  similarly  to  the  value  of  X  in  the  previous  section.  By 
substituting  the  coordinates  (?,,£,  of  a  known  point  on  the  reloading  curve  into  Eq.  3.55; 
and  solving  for  7,  which  is  also  take  as  the  tangent  modulus,  Ec,  one  obtains: 


CT,    +  K 

e.+L 


Eq.  3.  56 
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Fig.  3-19.  Typical  loading,  unloading  and  reloading  in  compression 

When  concrete  goes  into  tension,  if  it  goes  back  to  compression  it  will  reload 
only  when  the  gap  resulting  form  the  previous  cycle  is  closed.  The  process  can  be 
visualized  in  Fig.  3-20. 


Fig.  3-20.  Concrete  behavior  with  gap 
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Bilinear  Model 

In  another  study,  Agrawal  et  al.  (1965)  claimed  that  the  response  of  under- 
reinforced  concrete  beams  is  governed  by  the  steel.  It  was  therefore  proposed  to  use  a 
simplification  of  the  concrete  model  to  simplify  the  analysis  of  doubly  reinforced 
concrete  beams.  Idealized  elastic  plastic  curves  were  drawn  with  a  yield  stress  equal  to 
the  nominal  strength  fc  value  for  concrete,  and  an  elastic  modulus  equal  to  the  average 
stiffness  of  the  initial  portion  of  the  actual  stress-strain  curve  as  shown  in  Fig.  3-21. 
Also,  the  tensile  strength  of  concrete  is  neglected  in  this  model.  In  FLPIER  the  stress- 
strain  curve  for  concrete  is  based  on  13  points  as  shown  in  Fig.  3-22.  Referring  to  Fig. 
3-22,  the  value  of  the  average  elastic  modulus  is  taken  as 


b(5) 


Eq.  3.  57 


and  the  value  of  the  yielding  strain  ey  is  taken  as 
„    _<*(2) 


fc 


Eq.  3.  58 


Fig.  3-21 .  Bilinear  model  for  concrete 
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*0'  Concrete 


Fig.  3-22.  Stress-strain  curve  for  concrete  in  FLPIER 


Strain  Rate  Effect 


The  strength  of  reinforced  concrete  sections  is  very  dependent  on  the  strain  rate 
(Shkolnik  (1996)),  and  degree  of  confinement  of  the  section  (Saatcioglu  and  Razvi 
(1992),  and  Soroushian  et  al.  (1986)).  At  a  fast  strain  rate,  both  the  modulus  of  elasticity 
and  the  strength  of  concrete  increase.  The  increase  in  strength  can  be  as  much  as  17  % 
for  a  strain  rate  of  0.01/sec,  as  reported  by  Park  and  Paulay  (1975).  It  is  very  difficult  to 
test  structural  members  under  such  conditions.  So  before  using  the  results  from  static 
tests,  it  is  important  to  consider  the  effect  of  strain  rate  in  dynamic  analysis.  Important 
findings  about  the  strain  rate  effects  on  reinfoced  concrete  members  are  summarized 
next  (Otani,  (1980)): 
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1)  High  strain  rates  increased  the  initial  yield  resistance,  but  caused  small 

differences  in  either  stiffness  or  resistance  in  subsequent  cycles  at  the  same 
displacement  amplitude. 

2)  Strain  rate  effect  on  the  resistance  diminished  with  increased  deformation  in  a 

strain-hardening  range. 

3)  Non  substantial  changes  were  observed  in  ductility  and  overall  energy 

absorption  capacity. 
Otani  (1980)  also  suggests  that  the  strain  rate  during  an  oscillation  is  highest  at 
low  stress  levels,  and  gradually  decreases  toward  a  peak  strain.  Cracking,  crushing  and 
yielding  will  contribute  to  a  reduction  in  the  system's  stiffness,  elongating  the  period  of 
oscillation.  Such  damage  is  caused  by  the  lower  modes  of  vibration  having  long  periods. 
Therefore,  the  strain  rate  effect  is  small  in  earthquake  analysis,  and  has  a  small  effect  on 
the  response.  So,  the  static  hysteretic  behavior  can  be  used  in  a  nonlinear  dynamic 
analsyis  of  reinforced  concrete  structures. 

Confinement  F.ffect 

The  confinement  of  reinforced  concrete  columns  is  a  good  way  to  increase  its 

strength  and  ductility,  as  shown  in  Fig.  3-23. 

Confined  column 
\  

Moment 

Unconfined  column 


Curvatute 
Fig.  3-23.  Confined  and  unconfined  concrete  models  response 
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Confinement  basically  decreases  the  slope  of  the  descending  branch  of  the 
loading  curve  of  concrete,  making  the  confined  concrete  member  more  ductile  and  less 
brittle.  The  reader  is  referred  to  Saatcioglu  and  Razvi  (1992),  and  Soroushian  et  al. 
(1986)  for  a  good  discussion  on  the  subject. 

Soroushian  et  al.  (1986)  proposed  a  very  simplified  model,  that  incorporates 
both,  strain-rate  and  confinement  effects  in  the  compression  envelop  curve  for  concrete. 
The  following  constitutive  model  (Soroushian  et  al.  (1986))  was  implemented  in 
FLPIER: 


/  = 


*.*,/; 


2e 


^i^/c'fl  - z(e  - 0.002/:,^ )] 

e>  0.002  a:,  a:3 


where 


/=  concrete  compressive  stress 
£  =  concrete  compressive  strain 

P.f* 


K,=l  +  - 


./;" 


Eq.  3.  59 


4A„ 


ps  =  volumetric  ratio  of  the  hoop  reinforcement  to  concrete  core  =  — —  =  — 2. 
f'tm  28-day  compressive  strength  of  concrete,  adopted  as  0.85/'c 
fsh  =  yield  strength  of  transverse  reinforcement 
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29/; 
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(MPa) 


Eq.  3.  60 


h '  -  width  of  concrete  core  measured  to  outside  of  the  transverse  reinforcement, 
as  shown  in  Fig.  3.24  for  square,  rectangular,  and  circular  sections. 


\i            o 
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< ► 

h- 


■o                 O                O 
O                   n                 r. 

h1 

Fig.  3-24.  Core  width  for  different  cross  sections 


s  -  center-to  center  spacing  of  transverse  reinforcement 

K2=  1.48  +  0.160  log10E+0.0127(logl0e)2 
A:3=1.08  +  0.1601og,0E+0.0193(logl0e)2 


Eq.  3.  61 
Eq.  3.  62 


Note:  For  8  <  1 0"05  /  sec,  K2  =  K3  =  1 .0. 

In  this  model,  K2  represents  the  strain  rate  effect  on  the  compressive  strength  of 
the  concrete,  and  K}  takes  care  of  the  strain  rate  effect  on  the  strain  at  maximum  stress. 
It  is  assumed  that  the  strain  rate  effect  on  the  slope  of  the  descending  branch  of  the 
stress-strain  diagram  is  similar  to  the  strain  rate  effect  on  the  compressive  strength  of 
concrete.  This  is  supported  by  test  results.  The  model  representation  can  be  seen  in  Fig. 
3-25  below.  There  is  no  change  in  the  unloading  or  reloading  curves. 
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Concrete  Compressive  Stress  (0 


K,Kfc 


0.2K,Kf, 


0.002KiKi  Concrete  Compressive  Strain  (e) 

Fig.  3-25.  Confined  concrete  model 


The  following  modifications  were  also  proposed  by  Soroushian  et  al.  (1986)  for 
the  secant  and  tangent  stiffness  of  concrete  when  subjected  to  dynamic  loading: 


^-  =  1.241  +  0.1111og10E+0.127(log,0E)2 


-f-  =  1 .061  +  0.464  log10  e+  0.00683(logl0  s)2 


Eq.  3.  63 
Eq.  3.  64 


Where  Ecd  =  dynamic  secant  modulus  of  elasticity,  Ecs  =  static  secant  modulus  of 
elasticity,  E,d  =  dynamic  tangent  modulus  of  elasticity,  and  E,s  =  static  tangent  modulus 
of  elasticity.  Note  that  when  compared  to  the  secant  modulus,  the  tangent  modulus 
seems  to  be  less  influenced  by  the  rate  of  straining.  These  changes  were  also 
implemented  in  FLPIER. 


CHAPTER  4 
MODAL  ANALYSIS 

Because  the  most  frequently  used  procedure  for  designing  bridges  for 
earthquakes  is  based  on  modal  analysis,  a  brief  explanation  of  this  method  of  analysis  is 
given  next.  It  should  be  noted  that  modal  analysis  is  an  extensive  subject  and  that  what 
is  presented  here  is  only  an  introduction  of  the  basic  concepts  necessary  to  understand 
the  method.  The  interested  reader  is  referred  to  Chopra  (1995),  or  Paz  (1985),  for  a  more 
detailed  discussion  on  modal  analysis. 

Natural  Vibration  Frequencies  and  Modes 

Before  getting  to  modal  analysis  it  is  opportune  to  introduce  the  eigenvalue 
problem  whose  solution  gives  the  natural  frequencies  and  vibration  modes  of  a  system. 
The  free  vibration  of  an  undamped  system  in  one  of  its  natural  vibration  modes  can  be 
described  mathematically  by 

d(t)  =  q„m„  Eq.4.1 

where  the  deflected  shape  <)>„  does  not  vary  with  time.  The  time  variation  of  the 
displacements  is  described  by  the  simple  harmonic  function 

q„ (t)  =  An  cosco,/  +  Bn  sinconr  Eq.  4.  2 

where  A„  and  B„  are  constants  of  integration  that  can  be  determined  from  the  initial 
conditions  that  initiate  the  motion.  Combining  Eqs.  4.1  and  4.2  gives 
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d(0  -  •!>,,  (A  cosmj  +  Bn  sinco,/)  Eq.  4.  3 

where  co„  and  $p„  are  unknown. 

Substituting  this  form  of  d(t)  in  the  undamped  equation  of  motion  gives 

[-<°><P„  +*♦„]?„(')  =  0  Eq.4.4 

This  equation  can  be  satisfied  in  one  of  the  two  ways.  Either  q„(t)  =  0,  which 
implies  d(t)  =  0  and  there  is  no  motion  of  the  system  (this  is  the  so-called  trivial 
solution),  or  the  natural  frequencies  co„  and  modes  <|>„  must  satisfy  the  following 
algebraic  equation: 

**.=»»*♦»  Eq.  4.  5 

which  provides  a  useful  condition.  This  algebraic  problem  is  called  the  matrix 
eigenvalue  problem  (01  the  generalized  eigenvalue  problem).  The  stiffness  and  mass 
matrices  *  and  m  are  known;  the  problem  is  to  determine  the  scalar  co„2  and  the  vector 

To  indicate  the  formal  solution  to  Eq.  4.5,  it  is  rewritten  as 

[i -co  >]<|>„  =0  Eq.4.6 

which  can  be  interpreted  as  a  set  of  N  homogeneous  algebraic  equations  for  the  N 
elements  $M  (j=l,2,...^V).  This  set  always  has  the  trivial  solution  <|>„=0,  which  is  not 
useful  because  it  implies  no  motion.  It  has  nontrivial  solutions  if 

det\k-B>*m\  =  0  Eq.  4.  7 
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When  the  determinant  is  expanded,  a  polynomial  of  order  N  in  (B„2  is  obtained. 
Eq.  4.7  is  known  as  the  characteristic  equation  or  frequency  equation.  This  equation  has 
N  real  and  positive  roots  for  co„2  because  m  and  k,  the  structural  mass  and  stiffness 
matrices  are  symmetric  and  positive  definite  1.  The  positive  definite  property  of  k  is 
assured  for  all  the  structures  supported  in  a  way  that  prevents  rigid-body  motion.  Such  is 
the  case  for  civil  engineering  structures  of  interest  to  us,  but  not  for  unrestrained 
structures  such  as  aircraft  in  flight.  The  positive  definite  property  of  m  is  also  assured  if 
m  is  consistent,  and  can  be  assured  when  m  is  lumped  by  having  positive  masses  in  all 
diagonal  terms. 

The  N  roots  of  Eq.  4.7  determine  the  N  natural  frequencies  co„  (n  =  1 ,2,. .  .JV)  of 
vibration.  These  roots  of  the  characteristic  equation  are  also  known  as  eigenvalues, 
characteristic  values,  or  normal  values.  When  a  natural  frequency  <b„  is  known,  Eq.  4.6 
can  be  solved  for  the  correspondent  vector  <|)„  to  within  a  multiplicative  constant.  The 
eigenvalue  problem  does  not  fix  the  absolute  amplitude  of  the  vectors  <|>„,  only  the  shape 
of  the  vector  given  by  the  relative  values  of  the  N  displacements  fy„  (/'  =  1,2,...,JV). 
Corresponding  to  the  N  natural  vibration  frequencies  co„  of  an  jV-DOF  system,  there  are 
N  independent  vectors  <|>„  which  are  known  as  natural  modes  of  vibration,  or  mode 
shapes  of  vibration.  These  vectors  are  also  known  as  eigenvectors,  characteristic 
vectors,  or  normal  modes. 


1  In  mathematical  terms  a  matrix  A  e  R"""  is  positive  definite  if  x7 Ax  >  0  for  all 
nonzero  x  e  R" .  In  particular  all  the  diagonals  entries  of  A  are  positive. 
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Modal  and  Spectral  Matrices 

The  A7  eigenvalues,  N  natural  frequencies,  and  A' natural  modes  can  be  assembled 
compactly  into  matrices.  Let  the  natural  mode  <(>„  corresponding  to  the  natural  frequency 
a>„  have  elements  ty„,  where  j  indicates  the  DOF's.  The  A' eigenvectors  then  can  be 
displayed  in  a  single  square  matrix,  each  column  of  which  is  a  natural  mode: 


<*>  =  [♦>]  = 


4b 

4b 


Eq.  4.  8 


The  matrix  O  is  called  the  modal  matrix  for  the  eigenvalue  problem,  Eq  4.5.  The 
N  eigenvalues  co„  can  be  assembled  into  a  diagonal  matrix  Q2,  which  is  known  as  the 
spectral  matrix  of  the  eigenvalue  problem,  Eq.  4.5: 


Q2 


Eq.  4.  9 


Each  eigenvalue  and  eigenvector  satisfies  Eq.  4.5,  which  can  be  rewritten  as  the 
relation 

kifn  =  H»„a>„2  Eq  4_  10 

By  using  the  modal  and  spectral  matrices,  it  is  possible  to  assemble  all  of  these 
relations  into  a  single  matrix  equation: 


£®  =  wOO2 


Eq.4.  11 
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Equation    4. 1 1  provides  a  compact  presentation  of  the  equations  relating  all 
eigenvalues  and  eigenvectors. 

The  orthogonality  of  natural  modes  implies  that  the  following  square  matrices 
are  diagonal: 

K  =  <S>Tk<$  M  =  <S>'m<i>  Eq.  4.  12 

where  the  diagonal  elements  are 

&n  =  kn  *♦«  Mn  =  §*m§n  Eq.  4. 13 

Since  m  and  k  are  positive  definite,  the  diagonal  elements  of  K  and  M  are 
positive.  They  are  related  by 

^„=w„2M„  Eq.  4. 14 

Normalization  of  Modes 

As  mentioned  earlier,  the  eigenvalue  problem,  Eq.  4.5,  determines  the  natural 
modes  to  only  within  a  multiplicative  factor.  If  the  vector  <j>„  is  a  natural  mode,  any 
vector  proportional  to  <j>„  is  essentially  the  same  natural  mode  because  it  also  satisfies 
Eq.  4.5.  Scale  factors  are  sometimes  applied  to  natural  modes  to  standardize  their 
elements  associated  with  amplitudes  in  various  DOF's.  This  process  is  called 
normalization.  Sometimes  it  is  convenient  to  normalize  each  mode  so  that  its  largest 
element  is  the  unity.  Other  times  it  may  be  advantageous  to  normalize  each  mode  so  that 
the  element  corresponding  to  a  particular  DOF,  say  the  top  floor  of  a  multistory 
building,  is  unity.  In  theoretical  discussions  and  computer  programs  it  is  common  to 
normalize  modes  so  that  the  M„  have  unit  values.  In  this  case 


M,  =  £">$.  =  1  O'ffKD  =  /  Eq.  4.  15 

where  /  is  the  identity  matrix,  a  diagonal  matrix  with  unit  values  along  the  main 
diagonal.  Equation  4.15  states  that  the  natural  modes  are  not  only  orthogonal  but  are 
normalized  with  respect  to  m.  They  are  then  called  a  mass  orthonormal  set.  When  the 
modes  are  normalized  in  this  manner,  Eqs.  4.13  and  4.14  become 

K„  =  ♦/*♦»  =  <o,J K  mm*  K  =  <D7'M>  -  O2  Eq.  4. 16 

Modal  Equations 

The  equations  of  motion  for  a  linear  MDOF  were  derived  in  Chapter  2  and  are 
repeated  here  for  convenience  in  local  coordinates: 

mti  +  cii  +  ku  =  p(t)  Eq.  4.  17 

As  mentioned  earlier,  the  displacements  u  of  an  MDOF  system  can  be  expanded 
in  terms  of  modal  contributions.  Thus  the  dynamic  response  of  a  system  can  be 
expressed  as 

N 

"(0  =  £<l>,?,.(0  =  <M')  Eq.4.  18 

Using  this  equation,  the  coupled  equations  4. 1 7  in  «,(f)  can  be  transformed  to  a 
set  of  uncoupled  equations  with  modal  coordinates  q„(l)  as  the  unknowns.  Substituting 
Eq.4. 18  in  Eq.4. 17  gives 

H  N  H 

X™f<iU0  +  £4r<7,(0  +  X><Mv(/)  =  .p(/)  Eq.4.  19 

r=\  r=l  r=l 

Premultiplying  each  term  in  this  equation  by  (J)/  gives 
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"  »  N 

S*«**r4fr(0  +  2l*.r«!*A(0+2]*.r«*,ffr(0-*.r/K0  Eq.4.20 

!■•'  r-\  „\ 

Because  of  the  orthogonality  relations  of  Eq.  4.20,  all  terms  in  each  of  the 
summations  for  m  and  k  vanish,  except  the  r  =  n  term,  reducing  this  equation  to 

N 

(+>*„  )<7„  CO  +  Z  (♦/«♦, )?,  (')  +  (*/*<(•»  )q„  (I)  =  tTP(t)  Eq.  4.  21 


^„<7„  +  £  C,„qrt  +  K„q„t  =  Pn  (?)  Eq.  4.  22 


where 


K  =  f><(>„    C„r  =  tfcif,      K,  =  <j)„7 %      i»  (/)  =  <$,,',  p(t)  Eq.  4.  23 

Equation  4.22  exists  for  each  n  =  1  to  N  and  the  set  on  N  equations  can  be 
written  in  matrix  form 

Mq  +  Cq  +  Kq  =  P(t)  Eq.  4.  24 

where  M  is  a  diagonal  matrix  of  the  generalized  modal  masses  M„,  A-  is  a  diagonal 
matrix  of  the  generalized  modal  stiffnesses  K„,  C  is  a  nondiagonal  matrix  of  coefficients 
Cnr,  and  P(r)  is  a  column  vector  of  the  generalized  modal  forces  Pn{f).  These  jV  equations 
in  modal  coordinates  are  coupled  through  the  damping  terms  because  Eq.  4.22  contains 
more  than  one  modal  velocity. 

The  modal  equations  will  be  uncoupled  if  the  system  has  classical  damping 
(Rayleigh  is  a  form  of  classical  damping).  For  such  systems  Cm  =  0  if  n  *  r  and  Eq.  4.22 
reduces  to 
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K4.  +  C„q„  +  K„q„  =  P„(t)  Eq.  4. 25 

where  C„  is  given  by 

c„  =+/«♦«  Eq.4.26 

This  equation  governs  the  response  of  the  SDOF  system  shown  in  Fig.  4-1. 
Dividing  Eq.  4.25  by  M„  gives 

>.(0 


q„  +2£„(B„?„  +co„2<?„ 


Af. 


Eq.  4.  27 


where  t,„  is  the  damping  ratio  for  the  nth  mode.  The  damping  ratio  is  usually  estimated 
on  experimental  data  for  structures  similar  to  the  one  being  analyzed.  Equation  4.27 
governs  the  nth  modal  coordinate  q„(t),  and  the  parameters  Mn,  Kn,  C„  and  P„(t)  depend 
only  on  the  nth-mode  <j>„  ,  not  on  other  modes.  Thus  we  have  N  uncoupled  equations  like 
Eq.  4.25,  one  for  each  natural  mode.  In  summary,  the  set  on  N  coupled  differencial 
equations  4.17  in  nodal  displacements  »/')  has  been  transformed  to  the  set  of  JV 
uncoupled  equations  (Eqs.  4.25)  in  modal  coordinates  q„(t). 


I" 


"fr 


K„ 


1n(0 


H, 


Pn(t) 


^s\\\^\\^SS\\\\\\\^^^^\^J?- 


Fig.  4-1.  Generalized  SDF  system  for  the  nth  natural  mode 
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For  given  dynamic  forces  defined  by  p(t),  the  dynamic  response  of  a  MDOF 
system  can  then  be  determined  by  solving  Eq.  4.25  or  4.27  for  the  modal  coordinates 
q„(t).  Each  modal  equation  is  of  the  same  form  as  the  equation  of  motion  for  a  SDOF 
system.  Thus  the  solution  methods  and  results  available  for  SDOF  systems  can  be 
adapted  to  obtain  solutions  q„{t)  for  the  modal  equations.  Once  the  modal  equations 
have  been  determined,  Eq.  4.18  indicates  that  the  contribution  of  the  nth  mode  to  the 
displacement  u(t)  is 

"„(')  =  *„  (')<?„(<)  Eq.4.28 

and  combining  these  modal  contributions  gives  the  total  displacements: 

"(0  =  £«„(')  =  Z<tW„(0  Eq.4.29 

n-1  n-l 

Calculating  all  <(>„  in  Eq.  4.6  seems  to  be  a  prohibitive  computational  expense  for 
large  systems.  But  for  many  problems,  the  higher-frequency  modes  participate  little  in 
the  structural  response  and  therefore  only  a  small  number  of  low-frequency  modes  need 
to  be  used.  Thus  only  the  first  m  equations  of  Eq.  4.6,  where  typically  m  «  N,  are 
solved  and  the  transformation  is  approximated  by 

"0)  =  Z«„(0  =  !>„?„(')    where  m<N  Eq.4.30 

This  procedure  is  known  as  classical  modal  analysis  or  the  classical  mode 
superposition  method  because  individual  (uncoupled)  modal  equations  are  solved  to 
determine  the  modal  coordinates  qn{t)  and  the  modal  responses  u„{t),  and  the  latter  are 
combined  to  obtain  the  total  response  u(t).  More  precisely,  this  method  is  called  the 
classical  mode  displacement  superposition  method  because  modal  displacements  are 
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superposed.  For  brevity  we  usually  refer  to  this  procedure  as  modal  analysis.  This 
analysis  method  is  restricted  to  linear  systems  with  classical  damping.  Linearity  of  the 
system  is  implicit  in  using  the  principle  of  superposition,  Eq.  4.29.  Damping  must  be  of 
the  classical  form  in  order  to  obtain  modal  equations  that  are  uncoupled,  a  central 
feature  of  modal  analysis. 

Element  Forces 

The  element  forces  can  be  obtained  by  two  procedures.  In  the  first  procedure  the 
nth-mode  contribution  r„(t)  to  an  element  force  r(t)  is  determined  from  nodal 
displacements  «„(/)  using  element  stiffness  properties.  Then  the  element  force 
considering  contributions  of  all  modes  is 

r<t)  =  /Lr.(0  Eq.4.31 

In  the  second  procedure,  the  equivalent  static  forces  associated  with  the  nth- 
mode  response  are  defined  using  the  relationship /„(()  =  k„u„(t),  and  using  Eq.4.16  gives 

/„(/)=  m  ><t>„?„(0  Eq.4.32 

Static  analysis  of  the  structure  subjected  to  these  external  forces  at  each  time 
instant  gives  the  element  force  r„(t),  then  the  total  force  r(t)  is  given  by  Eq.  4.3 1 . 

Modal  Equations  for  Ground  Motion 

The  equations  of  motion  governing  the  response  of  an  MDF  system  to 
earthquake  induced  ground  motion  are  repeated: 

mu  +  cu  +  ku  =  pj(t)  Eq.4.33 
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where 

peB(t)  =  -miug(t)  Eq.4.34 

The  displacement  u  of  an  N-DOF  system  can  be  expressed  as 

N  N 

"C)  =  I>„W  =  I>„9„(0  Eq.4.35 

The  spatial  distribution  of  the  effective  earthquake  forces  pejj(t)  is  defined  by  s  = 
m\  .  this  force  distribution  can  be  expanded  as  a  summation  of  modal  inertia  force 
distributions  s„: 

N 

2XH>„  Eq.4.36 


OTl  =  ?_.l  „m 


where 


r»=M"  **-*>  M.=tfn*,  Eq.4.37 


the  factor  f„  is  usually  called  the  modal  participation  factor,  implying  it  is  a  measure  of 
the  degree  to  which  the  nth  mode  participates  in  the  response.  The  contribution  of  the 
nth  mode  to  the  excitation  vector  m\  is 

*»  =  r„H>„  Eq.  4.  38 

which  is  independent  of  how  the  modes  are  normalized. 

Equation  4.27  is  then  specialized  for  earthquake  motion  by  replacing  p(t)  by 
Peg(t)  to  obtain 

9.  +2£„w„?„  +o„2(sf„  =  -!*„«,  (0  Eq.  4.  39 
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The  solution  of  q„(t)  can  readily  be  obtained  by  comparing  Eq.  4.39  to  the 
equation  of  motion  for  the  nth-mode  SDF  system,  an  SDOF  system  with  vibration 
properties-natural  frequency  co„  and  damping  ratio  &,  -  of  the  nth  mode  of  the  MDOF 
system.  Equation  4.39  with  g  =  \„  ,  which  governs  the  motion  of  this  SDOF  system 
subjected  to  ground  motion  itg(t) ,  is  shown  here  with  u  replaced  by  D„  to  emphasize  its 
connection  to  the  nth  mode: 

A,  +2§„(D„D,,  +»,'£>,,  =  -iig(t)  Eq.  4.  40 

Comparing  Eq.  4.39  to  4.40  gives 

?„(r)  =  r„D„(0  Eq.4.41 

Thus  q„(t)  is  readily  available  once  Eq.  4.40  has  been  solved  for  D„(t),  utilizing 
numerical  time  stepping  methods  for  SDOF  systems. 

As  mentioned  earlier  there  are  two  procedures  available  to  determine  the  forces 
in  the  elements  from  the  displacements  u„(t).  The  second  of  these  procedures,  using 
equivalent  static  forces,  is  preferred  in  earthquake  analysis  because  it  facilitates 
comparison  of  dynamic  analysis  procedures  with  earthquake  design  forces  specified  in 
building  codes.  Equation  4.32  defines  the  static  forces  associated  with  the  nth-mode 
response,  where  q„(i)  is  given  by  Eq.  4.41.  Putting  these  equations  together  and  using  s„ 
(Eq.  4.38)  leads  to 

/„(')  =  •?„ A„(t)  Eq.4.  42 

where, 

A„(t)  =»,2Z).(0  Eq.4.  43 
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The  resultant  equivalent  static  forces  f„(t)  are  the  product  of  two  quantities:  (1) 
the  nth-mode  contribution  s„  to  the  spatial  distribution  mi  of  pej(t),  and  (2)  the  pseudo- 
acceleration  response  of  the  nth-mode  SDOF  system  to  u  (I) . 

Response  Spectrum  Analysis 

The  type  of  modal  analysis  that  has  been  presented  is  called  response  history 
analysis  (RHA),  because  the  response  depends  on  analyzing  N  uncoupled  SDOF 
equations  of  motion  on  time  by  any  procedure  used  for  SDOF  systems.  However 
structural  design  is  usually  based  on  the  peak  values  of  forces  and  deformations  over  the 
duration  of  the  earthquake-induced  response.  Therefore  another  type  of  modal  analysis 
was  developed,  which  is  called  response  spectrum  analysis  (RSA). 

This  method  is  based  on  the  fact  that  the  exact  value  of  the  peak  response  of  an 
MDOF  system  in  its  nth  natural  mode  can  be  obtained  from  the  earthquake  response 
spectrum.  A  spectrum  is  a  plot  of  the  period  T  or  natural  frequency  m  versus  the 
maximum  (or  peak)  response,  which  could  be  acceleration,  velocity  or  displacement  for 
a  SDOF  system  subjected  to  a  specific  earthquake  with  a  specific  damping  ratio.  Such  a 
plot  is  shown  in  Fig.  4-2.  Therefore  the  peak  value  of  A„(t)  is  available  from  the  pseudo- 
acceleration  response  spectrum  as  its  ordinate  A(T„£„)  corresponding  to  the  natural 
period  T„  and  damping  ration  £„;  for  brevity,  A(T„£„)  will  be  denoted  as  A„.  Therefore 
the  peak  response  is 


r»  =r'A. 


Eq.  4.  44 
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All  response  quantities  r„(t)  associated  with  a  particular  mode,  say  the  nth  mode, 
reach  their  peak  values  at  the  same  time  instant  as  A„(t)  reaches  its  peak. 


Maximum  u,  v 
or  a 


co  ot  T 
Fig.  4-2.  Typical  response  spectrum     "s  (8lven) 


Modal  Combination  Rules 

How  can  we  combine  the  peak  modal  responses  rno  (n  =  1,2,...^)  to  determine 
the  peak  value  r0s  max  |  r(t)  \  of  the  total  response?  It  will  not  be  possible  to  determine 
the  exact  value  of  r0  from  rno  because,  in  general,  the  modal  responses  r„(t)  attain  their 
peaks  at  different  time  instants  and  the  combined  response  r(f)  attains  its  peak  at  yet  a 
different  instant. 

Approximations  must  be  introduced  in  combining  the  peak  modal  responses  r„„ 
determine  from  the  earthquake  response  spectrum  because  no  information  is  available 
when  these  peak  modal  values  occur.  The  assumption  that  all  modal  peaks  occur  at  the 
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same  time  and  their  algebraic  sign  is  ignored  provides  an  upper  bound  to  the  peak  value 
of  the  total  response: 

/-„<XkJ  Eq.4. 45 

w=1 

This  upper-bound  value  is  usually  too  conservative.  Therefore,  this  absolute  sum 
(ABSSUM)  modal  combination  rule,  is  not  popular  in  structural  design  applications. 
The  square-root-of-sum-of-squares  (SRSS)  rule  for  modal  combination  is 


The  peak  response  in  each  mode  is  squared,  the  squared  modal  peaks  are 
summed,  and  the  square  root  of  the  sum  provides  an  estimate  of  the  peak  total  response. 
This  modal  combination  rule  provides  excellent  response  estimates  for  structures  with 
well-separated  frequencies. 

The  complete  quadratic  combination  (CQC)  rule  for  modal  combination  is 
applied  to  a  wider  class  of  structures  at  is  overcomes  the  limitations  of  the  SRSS  rule. 
According  to  the  CQC  rule, 

(UN  \  "2 

r. 


EEp-V-  Eq.4. 47 

Each  of  the  N1  terms  on  the  right  side  of  this  equation  is  the  product  of  the  peak 
response  in  the  ith  and  nth  modes  and  the  correlation  coefficient  pin  for  these  two 
modes;  p,„  varies  between  0  and  1  and  p,„  =  1  for  i  =  n.  Thus  Eq.  4.47  can  be  rewritten 
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N     N 

.:„        J 

Eq.  4.  48 


to  show  that  the  first  summation  on  the  right  side  is  identical  to  the  SRSS  combination 
rule  of  Eq.  4.46;  each  term  in  this  summation  is  obviously  positive.  The  cross-modal 
coefficient  p,„  can  be  approximated  by 

(l "  Pi )  +  «&.  P„,  (l  +  P^  )  +  4fe2  +  4,2  )P,' 

where  p/n  =  ai  I  ta„.  This  equation  also  implies  that  p,„  =  pnl,  p,„  =  1  for  ;'  =  n  or  for 
two  modes  with  equal  frequencies  and  equal  damping  ratios.  For  equal  modal  damping 
%t  =  'ti„  =  \  this  equation  simplifies  to 

842(i  +  P(JP,f 

The  topic  of  modal  analysis  has  been  presented  in  a  very  condensed  and 
simplified  manner.  The  subject  is  so  extensive  that  a  whole  volume  could  be  written  on 
it.  The  interested  reader  should  consult  a  structural  dynamics  book  (Chopra  (1995),  Paz 
(1985),  Paz  (1994),  Craig  (1981),  among  others)  for  a  more  detailed  discussion.  It  is 
always  important  to  emphasize  that  modal  analysis  relies  on  the  hypotheses  that  the 
superposition  of  the  responses  of  SDOF  systems  is  used  to  get  the  response  of  a  MDOF 
system.  This  implies  that  the  mass  and  stiffness  matrices  remain  constant  during  the 
process,  in  other  words,  a  linear  analysis  is  always  performed.  Furthermore  if  damping 
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is  considered,  the  damping  matrices  must  not  only  remain  constant  during  the  process, 
but  also  be  in  the  classical  form  (i.e.  proportional  to  the  mass  and  stiffness  matrices). 

How  FLPIER  handles  Modal  Analysis 

The  computer  program  FLPIER  has  also  been  implemented  with  modal  analysis 
capabilities,  referring  to  Fig.  4-3  the  following  iterative  procedure  was  adopted: 
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Fig.  4-3.  Modal  analysis  of  pier 


In  the  first  cycle  the  earthquake  is  applied  to  the  structure  and  the  initial  forces  at 
the  base  of  the  piers  are  computed.  Initially  the  springs  that  represent  the  foundation  are 
considered  very  stiff,  to  simulate  fixed  supports.  Then  for  each  column,  a  vector  of  six 
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forces  is  generated,  the  three  forces  Fx,  Fy,  and  Fz  in  the  x,  y,  and  z  directions,  and  the 
three  Mx,  My,  and  Mz  respective  moments. 

Then  each  of  these  forces  is  applied  to  the  foundation,  one  at  a  time,  like  in  a 
regular  static  analysis.  This  will  produce  three  displacements,  dx,  dy,  and  dz,  in  the  x,y, 
and  z  directions,  and  three  respective  rotations  6y,  0y,  and  9z,  at  the  base  of  each 
column.  These  define  the  first  column  of  the  flexibility  matrix  for  the  foundation.  After 
all  six  forces  are  applied  we  have  the  six  by  six  flexibility  matrix  for  the  foundation,  one 
for  each  column.  Inverting  this  matrix  we  obtain  the  new  stiffness  for  the  foundation, 
which  becomes  the  foundation  springs  for  the  base  of  each  pier  for  the  next  cycle. 

The  analysis  is  carried  out  until  two  consecutive  base  forces  for  one  pier  are  the 
same  within  a  tolerance.  To  compare  vectors  we  use  norms,  also  called  absolute  values. 
The  2-norm2  of  two  consecutive  base  force  vectors  is  computed,  and  if  within  a 
stipulated  tolerance,  the  analysis  is  terminated  and  the  final  forces  are  printed. 

It  is  important  to  note  that  in  this  analysis  the  structure  is  always  considered  to 
be  linear,  the  main  pre-requisite  for  modal  analysis  of  any  type.  However,  the  springs 
generated  for  each  cycle  will  reflect  the  characteristic  nonlinear  behavior  of  the 
foundation.    This    is    an    approximate    method    used    for    bridge    pier    analsyis. 


2  The  2-norm  of  a  vector  is  defined  as:  If  5  =  (ul...u„ )  then  the  2-norm  of  u  is 


CHAPTER  5 
MULTIPLE  SUPPORT  EXCITATION 

Usually  we  assume  that  all  supports  where  the  structure  is  connected  to  the 
ground  undergo  identical  prescribed  motion.  In  this  section  we  generalize  the  previous 
formulation  of  the  equation  of  motion  to  allow  different  -  possibly  even  multi- 
component  -  prescribed  motions  at  the  various  supports.  Such  multiple-support 
excitation  may  arise  in  several  situations.  First,  consider  the  earthquake  analysis  of 
extended  structures  such  as  the  Golden  Gate  Bridge,  in  San  Francisco,  California.  The 
ground  motion  generated  by  an  earthquake  on  the  nearby  San  Andreas  fault  is  expected 
to  vary  significantly  over  the  6450-ft  length  of  the  structure.  Therefore  different  motions 
should  be  prescribed  at  the  four  supports:  the  base  of  the  two  towers,  and  two  ends  of 
the  bridge.  Second,  consider  the  dynamic  analysis  of  bridge  foundations.  The 
earthquake  occurs  deep  in  the  soil,  at  the  rock  level.  As  the  shock  wave  moves  up  in  the 
soil  the  acceleration  records  at  different  soil  depths  are  expected  to  vary.  Therefore  the 
piles  will  be  subjected  to  different  accelerations  along  its  depth  during  the  earthquake. 

In  order  to  analyze  multiple  support  excitation  systems,  we  must  include  in  the 
equation  of  motion  the  effects  of  the  support  motions  (Fig.  5-1).  We  first  divide  the 
displacement  vector  into  two  parts  El  and  Dg  ;  £>'  includes  the  NDOF  of  the 
superstructure,  where  the  superscript  I  denotes  that  these  are  total  displacements;  and 
Dg  contains  the  Ng  components  of  support  displacements. 
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The  equation  of  dynamic  equilibrium  for  all  the  DOF  can  now  be  rewritten  in 
partitioned  matrix  form: 


D' 


c       cs 

el    c„„ 


*     ** 

kl    k„ 


0 


Eq.  5. 1 
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Superstructure 
DOF:^ 


Support  DOF:  D, 


Fig.  5.1.  Multiple  support  motion 


Observe  that  no  external  forces  are  applied  along  the  superstructure  DOF.  In  Eq. 
5.1  the  mass,  damping,  and  stiffness  matrices  can  be  determined  from  the  properties  of 
the  structure  using  the  traditional  methods  of  matrix  analysis,  while  the  support  motions 
°s(')>  ^»(0  and  ^g(0  must  be  specified.  Because  the  data  that  is  usually  available 
for  an  earthquake  is  the  acceleration  record  Dg(t),  the  quantities Dg(t)  and  Dg(t)  can 
be  obtained  by  numerical  integration  of  Dg(t),  by  using,  for  example,  the  Trapezoidal's 
rule.  It  is  desired  to  determine  the  displacements  D'  in  the  superstructure  DOF  and  the 
support  forces  Pg(t). 

Let's  first  separate  the  displacements  into  two  parts: 


D'\      \D" 
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D.I      ID  |+    0  E"-5-2 


In  this  equation  Lf  is  the  vector  of  structural  displacements  due  to  static 
application  of  the  prescribed  support  displacements  Dg  at  each  time  instant.  The  two  are 
related  through 


*    K 

kl    k,. 


D.rw 


where  p"g  are  the  support  forces  necessary  to  statically  impose  displacements  Ds  that 
vary  with  time;  obviously,  If  varies  with  time  and  therefore  is  known  as  the  vector  of 
quasi-static  displacements.  Observe  that  p'g  =  0  if  the  structure  is  statically  determinate 
or  if  the  support  system  goes  rigid-body  motion;  for  the  latter  condition  an  obvious 
example  is  identical  horizontal  motion  of  all  supports.  The  remainder  D  of  the  structural 
displacements  are  known  as  dynamic  displacements  because  a  dynamic  analysis  is 
necessary  to  evaluate  them. 

With  the  total  structural  displacements  split  into  quasi-static  and  dynamic 
displacements,  Eq  5.2,  we  return  to  the  first  of  the  two  partitioned  Eq.  5.1. 

mb'  +  mgDg  +  cD'  +cgDg  +  kD'  +kgDg=0  Eq.  5.  4 

Substituting  Eq.  5.2  and  transferring  all  terms  involving  Dg  and  Cf  to  the  right 
side  leads  to 

mb  +  cD  +  kD  =  peI(t)  Eq.  5.5 
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where  the  vector  of  effective  earthquake  forces  is 

p„(0  =  -(mD<  +mgDg)-(cD<  +cgDg)-(kD'  +  kgDg)  Eq.  5.  6 

This  effective  force  vector  can  be  rewritten  in  a  more  useful  form.  The  last  term 
drops  out  because  of  Eq.  5.3  resulting 

kD'+kgDg=0  Eq.5. 7 

This  relation  also  enables  us  to  express  the  quasi-static  displacements  If  in 
terms  of  the  specified  support  motions  Dg,  and  we  obtain 

D'  =  (Dg  Eq.  5.  8 

and 

l—k-\  Eq.5. 9 

We  call  I  the  influence  matrix  because  it  describes  the  influence  of  support 
displacements  on  the  structural  displacements.  If  all  the  supports  undergo  the  same 
motion  the  influence  matrix  I  is  unitary,  and  we  obtain  the  equations  derived  in  Chapter 
2.  It  is  clear  form  Eq,  5.9  that  t  is  obtained  by  solving  the  linear  system  of  equations: 

kt  =  -*,  Eq.  5.  10 

Substituting  Eqs.  5.8  and  5.7  in  Eq.  5.6  gives 

peff(t)  =  -(me  +  mg)De(l)-(ce  +  cg)Dg(t)  Eq.5. 11 
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If  the  ground  (or  supports)  accelerations  Dg(f)  are  prescribed,  the  velocities 
De(t)  can  be  obtained  by  numerical  integration  of  the  accelerations  Dg(t),  andp^t)  is 
known  from  Eq.  5.1 1,  and  this  completes  the  formulation  of  the  governing  equation  (Eq. 
5.5). 

Equation  5.1 1  can  be  simplfied  on  two  counts.  First,  if  the  damping  matrices  are 
proportional  to  the  stiffness  matrices  (i.e.  c  =  a/k  and  cg  =  aikg),  the  damping  term  is 
zero  because  of  Eq.  5.7.  Second,  if  the  mass  is  idealized  as  lumped  at  the  DOF,  the 
mass  matrix  is  diagonal,  implying  that  mg  is  a  null  matrix  and  m  is  diagonal.  With  these 
simplifications  Eq.  5.1 1  can  reduces  to 


ptI{t)  =  -mWg(l) 


Eq.  5. 12 


For  a  better  understanding  of  how  the  influence  t  matrix  is  formed,  consider  the 
2D  frame  with  the  DOF  illustrated  in  Fig.  5-2  below 

3 
1 


r-^—i 


5 


Fig.  5-2.  2D  frame  submitted  to  multiple  support  motion 
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Note  that  the  support  DOF,  4  and  5,  are  numbered  last,  according  to  Eqs.  5.1.  By 
applying  unit  displacements  to  DOF  4  and  5  we  get  the  4th  and  5th  columns  of  k  shown 
below 


EI 


24  6i  6L  12  12 

6L  U2  2L1  6L  0 

61  2L2  U2  0  61 

12  61  0  12  0 

12   0  61  0  12 


Eq.  5. 13 


If  we  make  kg  =  \kgl,kg2   ,  where  kgj  and  kgj  are  respectively  the  4th  and  5th 

columns  of  k,  and  the  influence  matrix  (  =  \( , ,  £  2  ]  (note  that  (  is  N  x  Ng),  by  solving 
H,  =  kgJlwe  get  «f/t  column  of  £.  This  procedure  can  be  extended  to  3D  analysis  as 
well.  In  this  case  there  are  six  possible  directions  for  the  motion  of  each  support, 
however  in  most  cases  the  horizontal  movement  controls  the  analysis.  In  FLPIER  the 
support  motions  are  applied  to  the  soil  springs.  In  the  case  of  the  simplified  pile  shown 
in  Fig.  5-3,  considering  the  DOF  illustrated,  the  support  stiffness  kg  (8  x  4)  is  then  given 


by  Eq.  5-14. 
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Fig.  5-3.  Pile  subjected  to  multiple  support  excitation 
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Eq.  5. 14 


Modal  analysis  has  also  been  developed  for  multi-support  excitation.  The  reader 
should  refer  to  Monti  et  al.  (1996),  Der  Kiurieghian  (1992),  Der  Kiurieghian  (1995), 
and  Chopra  (1995)  for  more  details. 


CHAPTER  6 
SOIL  STRUCTURE  INTERACTION 

The  earthquake  response  of  bridge  structures  depends  on  the  soil-pile-structure 
interaction  during  the  earthquake  loading.  The  subject  has  been  extensively  studied  by 
many  researchers  (Nogami  (1997),  Prakash  (1997),  Turner  (1995),  Bandoni  (1996), 
Reese  (1974),  Matlock  (1978),  Anderson  (1972),  Ting  (1987),  and  McVay  et  al. 
(1996)).  There  are  two  basic  limitations  in  these  studies.  The  first  one  is  the  fact  that 
most  of  them  were  done  for  single  pile  structures,  making  it  impossible  to  account  for 
multiple  pile  interaction  effect,  characteristic  of  a  real  structure.  The  second  one  is  that 
because  of  obvious  cost  problems,  large  scale  tests  are  rarely  done.  Another  limitation  is 
the  fact  that  the  original  data  for  most  of  these  studies  is  not  available. 

In  FLPIER,  the  soil-structure  interaction  can  be  accounted  for  by  the  use  of 
equivalent  nonlinear  springs.  In  particular,  the  complicated  soil-pile-superstructure 
interaction  is  captured  by  determining  the  primary  structural  earthquake  response,  and 
then  driving  this  response  into  the  foundation  system.  Thus,  the  coupled  problem  of  soil 
structure  interaction  is  solved  in  terms  of  structure  first  and  foundation  system  second, 
as  mentioned  in  chapter  5.  This  procedure  is  considerably  simpler  than  performing  a 
completely  coupled  solution  process;  however,  there  are  some  approximations  and 
uncertainty  in  the  simplified  procedure. 

Abghari  and  Chai  (1995)  suggest  that  the  main  source  of  uncertainty  in  the 
present  approach  can  be  illustrated  by  a  simple  modal  response  argument.  The  primary 
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response  of  the  structure  is  often  a  purely  structural  "fundamental"  mode  that  does  not 
capture  the  accurate  response  of  the  foundation.  The  primary  soil  response  mode 
approximates  a  free-field  shear  deformation  in  which  the  superstructure  is  displaced  in  a 
near-rigid-body  manner.  The  two  modes  may  not  be  closely  coupled  and  therefore 
driving  the  design  of  the  piles  with  the  response  of  the  superstructure  may  not  accurately 
estimate  the  actual  response  of  the  entire  system.  In  the  absence  of  more  rational 
analysis  data,  the  uncertainties  present  in  the  simplified  analysis  make  it  difficult  to 
predict  whether  the  resulting  pile  foundations  are  over-  or  under-designed.  It  is  therefore 
necessary  to  develop  models  that  can  account  for  the  coupling  between  foundation  and 
structure. 

A  simplified  coupled  solution  is  to  select  a  single  pile  from  the  pile  group  with 
the  superstructure  being  modeled  as  single  degree  of  freedom  system  having  the  same 
period  as  the  fundamental  period  of  the  structure,  and  a  mass  equal  to  the  contribution  of 
the  superstructure  to  each  pile.  The  pile  is  then  subjected  to  displacement  time  histories 
previously  calculated  from  a  dynamic  site  response  analysis  and  forces  and  defections 
are  calculated.  There  are  basically  two  methods  for  applying  the  displacement  histories 
to  the  piles:  the  Uncoupled  Method  and  the  Coupled  Method.  Both  methods  are 
described  next. 

Uncoupled  Method 

This  method  is  illustrated  in  Fig.  6-1  (a).  The  analysis  is  done  in  two  steps: 

1 )  The  computation  of  the  free  field  motions; 

2)  The  prediction  of  the  pile-superstructure  response. 
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The  free  field  motions  are  computed  independently  through  a  one-dimensional 
site  response  analysis  using  widely  available  computer  codes  (SHAKE.  Schnabel  et  al., 
1 972).  These  computed  ground  motions  are  then  used  as  boundary  conditions  applied  at 
nodal  locations,  corresponding  to  the  nonlinear  "p-y  springs",  in  the  foundation  system. 
Recalling  Chapter  5,  it  is  clear  that  what  is  done  is  a  multiple  support  excitation 
analysis,  with  the  support  motions  given  by  step  (1).  Usually  in  this  method  the  soil 
mass  is  considered  to  be  lumped  at  the  pile  nodes.  This  method  is  the  state-of-practice  in 
dynamic  analysis  of  foundations. 

Coupled  Method 

This  method  is  illustrated  in  Fig.  6-1  (b).  The  analysis  is  done  by  applying  the 
free  field  motion  to  the  base  of  a  soil  column,  modeled  with  "soil  elements".  The  idea  is 
that  the  earthquake  occurs  at  the  rock  level,  so  applying  the  free  field  motion  to  the 
bottom  of  the  soil  column  makes  sense.  The  "soil  elements"  are  connected  to  the  pile 
nodes  through  "p-y  springs"  and  dashpots.  Note  that  in  this  method  only  one  free  field 
motion  is  applied  to  the  structure-foundation  system,  so  there  is  no  necessity  for 
multiple  support  excitation.  One  advantage  of  this  method  is  that  the  soil  mass  can  be 
obtained  from  the  "soil  elements"  using  the  Finite  Element  Method  procedure  (e.g.  a 
consistent  mass  formulation).  The  disadvantage  of  this  model  is  that  the  "soil  element" 
must  be  a  good  representation  of  the  soil,  because  all  the  interaction  behavior  between 
the  support  motion  and  the  structure-foundation  system  is  dependent  on  this  element. 

Although  this  simplified  procedure  is  efficient  for  a  first  estimate  of  the  behavior 
of  the  structure  when  subjected  to  earthquake,  it  is  somewhat  impractical  for  design. 
Note  that  the  pile  cap,  pier  cap,  and  pile  group  effects  are  not  present  in  this  type  of 
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analysis,  making  it  difficult  to  estimate  the  forces  to  design  for.  Therefore  the  computer 
program  FLPIER  was  modified  to  perform  a  time-step  analysis  of  the  structure  using  the 
uncoupled  method.  It  is  now  possible  to  input,  for  each  pile  node,  the  acceleration 
record,  and  the  program  will  run  the  time-step  analysis  considering  a  multiple  support 
excitation.  All  piles  in  the  group  are  considered  to  be  submitted  to  the  same  acceleration 
record.  The  result  is  the  maximum  forces  that  each  element  has  been  subjected  to  during 
the  dynamic  analysis. 


Superstructure 


Superstructure 


U,i(t) 


U„<t) 


U„(t> 


(a)  (b) 

Fig.  6-1.  a)Coupled  model;  b)Uncoupled  model 
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Cyclic  Behavior  of  Soil 

The  cyclic  behavior  for  soil  presented  here  is  based  on  the  recommendations  of 
O'Neill  et  al.  (1997),  and  is  a  simplified  approach  to  the  complex  dynamic  behavior  of 
soil.  Although  the  changes  in  FLPIER  were  limited  to  the  lateral  resistance  springs  (p- 
y),  the  extension  to  the  axial  resistance  springs  (t-z)  is  immediate.  It  is  also  important  to 
note  that  the  dynamic  factor  usually  incorporated  in  the  p-y  curves  for  soil  under 
dynamic  loading  has  not  been  incorporated  in  this  model.  The  static  p-y  curves  are  used 
for  the  dynamic  analysis.  A  typical  p-y  curve  is  illustrated  in  Fig.  6-2. 

p(KN/m) 


y(m) 
Fig.  6-2.  Typical  p-y  curve 

The  p-y  behavior  is  illustrated  in  Fig.  6-3.  Loading  is  described  by  the  original  p- 
y  curve.  As  the  pile  loads  is  reversed  at  a  given  level,  the  soil  reaction  unloads  along  a 
path  parallel  to  the  loading  path,  not  along  the  initial  loading  path,  which  produces 
hysteretic  damping.  When  the  soil  reaction  reaches  zero,  no  further  reaction  is  generated 
until  the  displacement  of  the  pile  at  the  location  of  concern  reverses,  at  which  time 
follows  the  loading  path  in  the  opposite  direction.  Then,  upon  load  reversal  in  the 
opposite  direction,  a  mirror  image  effect  is  generated,  except  that  no  further  soil  reaction 
will  be  generated  until  the  defection  reaches  the  value  of  deflection  corresponding  to  the 
width  of  the  gap  in  the  previous  cycle.  Then,  loading  occurs  along  the  previous 
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unloading  path.  Once  the  gap  develops  at  a  particular  level  due  to  lateral  loading,  the 
axial  performance  of  the  pile  is  affected  because  the  development  of  shearing  resistance 
is  no  longer  permitted  at  that  location,  however  this  effect  is  not  taken  into  account  in 
the  actual  dynamic  version  of  FLPIER  at  this  time. 


p(KN/m) 


y(m) 


Fig.  6-3.  Cyclic  soil  model 


Cyclic  Dearadation 

Cyclic  degradation  is  specified  through  the  use  of  Eq.  6.1  (O'Neill  et  al.,  1997): 
Pc=(l-)-ipp-P*)+Pd  Eq.6. 1 

In  this  equation  X  is  a  soil  degradation  parameter  specified  by  the  user,  pc  is  the 
soil  resistance  (e.  g.,p  forp->>  curve  corresponding  to  a  given  value  of  deflection  y)  for 
the  current  cycle  of  loading,  pp  is  the  value  of  resistance  corresponding  to  the  present 
vale  of  y  on  the  previous  cycle  of  loading,  and  pd  is  the  fully  degraded  value  of  p  at  the 
present  value  of  y.  The  analysts  inputs  the  back  bone,  or  first  cycle,  p-y  curve  and  also  a 
fully  degraded /?->■  curve.  In  FLPIER  the  initial  p-y  curve  can  be  automatically  generated 
by  the  program  based  on  the  soil  properties.  The  fully  degraded  p-y  curve  is  given  by  a 
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degradation  factor,  given  by  the  user,  multiplied  by  the  initial  p-y  curve.  At  this  time  the 
user  can  not  supply  the  fully  degraded  p-y  curve.  In  a  seismic  event  or  for  an  extreme 
event  that  involves  impact  loading,  some  cyclic  degradation  may  occur,  but  the  fully 
degraded  value  may  never  be  reached.  A  typical  value  for  the  soil  degradation 
parameter,  X,  is  given  in  Eq.  6-2. 

(Jail 
A.  =  1-10^'""'"^  Eq.6.2 

Where  N$o  is  the  number  of  cycles  that  would  be  necessary  to  degrade  the  soil 

by  50  percent.  This  would  be  a  site-specific  parameter  that  wold  have  to  be  obtained  by 

appropriate  laboratory  tests  or  from  cyclic  lateral  pile  testing. 

Strain  Rate  Effect 

Usually  when  soil  is  subjected  to  dynamic  loading,  cyclic  degradation  occurs 
simultaneously  with  an  increase  in  the  apparent  soil  (axial  or  lateral)  resistance  caused 
by  rapid  rates  of  loading  that  occur  in  the  most  extreme  events.  Once  the  cyclic 
degradation  has  been  accounted  for,  the  rate  effect  is  computed  from  (O'Neill  et  al., 
1997): 

^-  =  l  +  F2log(^|  Eq.6.3 

where  p,  is  the  instantaneous  soil  resistance  considering  both  cyclic  and  loading  rate 
effects,  pc  is  the  resistance  considering  only  cyclic  loading,  tr  is  the  actual  rate  of  loading 
in  Hz  or  unit  of  distance  per  second,  /s  is  the  corresponding  rate  of  loading  appropriate 
for  standard  slow  cyclic  loading  (typically  0.01  to  0.1  Hz)  and  F2  is  a  soil  factor,  which 
can  be  taken  as  0.01  -  0.03  for  sand,  0.02  -  0.07  for  silts,  0.02  -  0.12  for  clays,  and  0.01  - 
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0.03  for  calcareous  soils.  Both  the  rate  effect  and  cyclic  degradation  models  are  built 
into  FLPIER. 

Radiation  Damping 

Radiation  damping  is  modeled  through  the  use  of  dashpots  having  constants  C 
(O'Neill  et  al.,  1997)  attached  to  the  pile  nodes,  where 

Ch=2D^-(v,+v)  Eq.6.4 

g 

where  C/,  is  for  horizontal  (p-y)  resistance  in  units  of  FT/L2,  D  is  the  pile  diameter,  y  is 

the  unit  weight  of  the  soil,  g  is  the  acceleration  of  gravity,  v3  is  the  shear  wave  velocity 

of  the  soil,  which  would  need  to  be  estimated  at  a  given  site,  and  v  is  a  velocity  in 

between  the  shear  and  compression  wave  velocities  of  the  soil.  A  lower  bound  solution 

for  Ca,  which  is  used  in  FLPIER,  is  given  by  taking  v  =  vs: 

Ch=4Dy-v,  Eq.6.5 

g 

A  similar  expression  can  be  developed  later  for  axial  loading. 


CHAPTER  7 
PREDICTIONS  OF  RESPONSE 

This  Chapter  provides  a  comparison  between  the  response  of  test  columns 
reported  in  the  literature  (Baron  and  Venkatesan  (1969),  Chen  and  Atsuta  (1973),  Hays 
(1975),  Chai  et  al.  (1991),  Bousias  et  al.  (1995),  and  Wilson  et  al.  (1997))  and  the 
theoretical  predictions  of  FLPIER.  The  tests  include  different  cross  sections,  such  as 
steel  sections,  circular  and  square  reinforced  concrete  sections.  When  the  data  is 
available,  the  monotonic  and  the  cyclic  tests  are  compared.  Three  examples  considering 
soil  are  also  presented,  but  the  structure  is  considered  to  remain  linear  during  the 
analysis.  The  results  for  the  various  tests  are  presented  next. 

Example  1  -  Steel  Section  1 

Although  this  work  is  mainly  about  the  behavior  of  reinforced  concrete 
structures  it  is  opportune  to  verify  the  behavior  of  a  steel  section  to  validate  the  adopted 
bilinear  steel  model.  AW  14  x  176  steel  section  that  was  first  studied  by  Baron  and 
Venkatesan  (1969),  and  later  used  by  Chen  and  Atsuta  (1973)  in  similar  studies,  is  used 
for  comparison.  The  steel  is  ASTM-A36,  with  yielding  stress  .ft  =  36  ksi,  and  Young's 
Modulus  Es  =  29000  ksi.  The  history  of  deformations  is  given.  No  axial  load  is  applied. 
The  W  section  is  modeled  as  a  single  steel  H-pile,  composed  of  16  nonlinear  discrete 
elements,  the  default  in  FLPIER,  and  is  fixed  at  the  base.  A  mass,  a  damper,  and  a 
spring  were  attached  to  the  top  of  the  pile  to  minimize  the  dynamic  effects,  resulting  in  a 
'pseudo-dynamic'  model  with  imposed  displacements.  The  FLPIER  model  and  the 
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dynamic  parameters  are  shown  in  Fig.  7-1 .  A  comparison  between  the  actual  test  and  the 
FLPIER  model  is  shown  in  Fig.  7-2.  Based  on  the  fact  that  just  a  few  points  are  reported 
in  the  original  work,  the  comparison  is  very  good. 


m  =  1  kip  s2/in 


rH 


k  =  10000  kip/ in 

1 


c  =  200  kip  s/in 


Fig.  7-1.  Example  1  computer  model 


-1  0  1 

Curvature  xlOOO(rads) 


Fig.  7-2.  Comparison  FLPIER  x  reference  for  cyclic  loading 


118 
Example  2  -  Steel  Section  2 

In  Example  2  the  geometric  nonlinearities  already  incorporated  in  FLPIER  are 
included  in  the  analysis,  with  the  presence  of  an  axial  load.  In  this  test  the  tangent 
stiffness  approach  added  to  the  code  is  compared  to  the  current  approach,  the  secant 
stiffness,  considering  the  p-A  effects.  The  computer  model  and  the  dimensions  for  the 
cross  section  are  illustrated  in  Fig.  7-3.  The  FLPIER  results  are  compared  to  the  results 
using  a  computer  program  that  also  considers  material  and  geometric  nonlinearities 
developed  by  Hays  (1975).  A  summary  of  all  the  results  for  this  test  can  be  seen  in 
Table  7-1.  The  comparison  between  FLPIER  and  Hays  (1975)  is  shown  in  Fig.  7-4  and 
7-5. 


Units  are  inches. 
fy  =  50  ksi,  E,  -  29000  ksi 


P  =  537.5  kips 


k=\x  10lu  kip/in 


^AA/Vf^ 


L- 15' -ISO" 


^^ 


Fig.  7-3.  Example  2  computer  model 
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Table  7-1. 

Results  for 

Example  2 

FLPIER 

Hays 
(tangent) 

Ax(in) 

Secant 

Tangent 

S(kips) 

M(kip.in) 

S(kips) 

M(kip.in) 

S(kips) 

M(kip.in) 

0.5 

3.5 

906.4 

3.5 

906.4 

3.5 

907 

1.0 

7.1 

1812.7 

7.1 

1812.7 

7.1 

1810 

1.5 

11.0 

2721.6 

10.8 

2717.8 

10.6 

2720 

2.0 

NC 

NC 

11.4 

3145.5 

11.5 

3140 

2.5 

10.9 

3297.8 

10.9 

3290 

3.0 

9.7 

3363.2 

9.7 

3360 

3.5 

8.4 

3393.6 

8.4 

3390 

4.0 

7.0 

3407.3 

7.0 

3400 

NC  =  No  convergence  achieved  for  the  analysis. 
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1.5  2  2.5  3 

Lateral  top  displacement  (in) 

Fig.  7-4.  FLPIER  x  Hays,  shear  force  comparison 


1.5  2  2.5  3 

Lateral  top  displacement  (in) 
Fig.  7-5.  FLPIER  x  Hays,  bending  moment  comparison 


The  following  observations  can  be  made  from  this  example: 

a)  As  expected  the  tangent  and  the  secant  approaches  give  the  same  response 
until  a  point  where  the  secant  approach  does  not  converge  (between  1.5  and  2.0  in). 
Beyond  this  point  FLPIER  is  in  good  agreement  with  the  results  given  by  Hays  (1975). 
This  indicates  that  the  new  tangent  approach  introduced  in  the  FLPIER  model  seems  to 
work,  and  in  this  particular  example  is  more  stable  than  the  secant  approach. 

b)  Second,  the  p-A  (or  second  order)  moment  effect  is  also  modeled,  and  is  clear 
in  Fig.  7-4.  Note  that  the  shear  force  starts  to  decrease  as  the  second  order  moment 
dominates  the  response.  Although  this  is  not  a  dynamic  test,  its  importance  relies  on  the 
fact  that  it  tests  the  tangent  stiffness  approach  used  in  the  dynamic  analysis. 
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Example  3  -  Circular  Reinforced  Concrete  Column  1 

The  third  example  analyzed  by  FLPIER  was  the  full-scale  flexure  column 
presented  in  Chai  et  al.  (1991).  The  reference  mentions  that  this  column  represented  the 
current  (1991)  ductile  design  for  bridge  columns.  Table  7-2  summarizes  the  parameters 
for  the  test  column  which  was  subjected  to  a  constant  axial  compression  force  of  1000 
kips  and  a  lateral  cyclic  displacement  of  increasing  amplitudes  until  failure  of  the 
column.  Unfortunately  the  data  for  the  cyclic  test  was  not  available  for  comparison.  The 
computer  model  is  identical  to  the  one  shown  in  Fig.  7-3,  except  for  the  cross  section 
and  length  L,  now  30  ft.  The  different  parameters  for  the  tests  are  shown  in  Table  7-3. 
All  tests  have  the  same  material  properties. 

Table  7-2.  Design  details  for  Example  3 


Diameter 

60" 

Height 

30' 

Cover  to  main  bar 

4" 

Concrete  strength  fc 
Modulus  of  elasticity  Ec 

5.2  ksi 
41 10  ksi 

Longitudinal  steel 

Yield  strength  ,/y 

Modulus  of  elasticity  E,  (adopted) 

25#14 
68.9  ksi 
29000  ksi 

Transverse  steel 
Yield  strength^ 

#5  spiral  at  3.5" 
71.5  ksi 

Axial  force 

1000  kips 

Table  7-3.  Mode 

parameters  for  Example  3 

File 

Stiffness 

Analysis  type 

Confinement 

Convergence 

circ2.in 

Secant 

Static 

No 

OK 

circ2 1  .in 

Tangent 

Static 

No 

OK 

circ22.in 

Tangent 

Dynamic 

No 

NO 

circ23.in 

Tangent 

Dynamic 

Yes 

OK 
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The  first  step  was  to  perform  an  incremental  static  analysis  using  the  secant  and 
the  tangent  methods.  The  static  test  uses  the  actual  uniaxial  material  model  presented  in 
Chapter  3.  The  dynamic  test  uses  the  proposed  concrete  model  also  presented  in  Chapter 
3.  The  comparison  between  FLPIER  and  the  test  is  shown  in  Fig.  7-6  and  7-7  for  the 
shears  and  moments  respectively.  From  Fig.  7.6  to  7.7  note  that  although  the  column 
strength  values  are  very  close,  the  slope  of  the  curve  is  slightly  different.  Three  factors 
may  have  contributed  to  that: 

a)  Because  the  original  work  did  not  provide  the  values  for  the  modulus  of 
elasticity  for  steel  and  concrete,  they  had  to  be  adopted.  For  the  concrete  the  empirical 
value  given  by  the  formula  (Meyer,  1996): 

Ec  =  57000777  Eq.  7. 1 

was  adopted.  For  the  steel  29,000  ksi  was  adopted  as  the  modulus  of  elasticity ,  which  is 
typical  for  steel. 

b)  The  use  of  a  discrete  curve,  defined  by  straight  lines  segments,  instead  of  a 
continuous  curve  (parabola)  for  the  concrete  model  may  also  have  introduced  some 
error. 

c)  A  phenomena  called  anchorage-slip,  which  is  explained  in  more  details  later, 
may  be  present.  Basically  this  phenomena  causes  additional  displacements  usually  not 
accounted  for  by  the  analysis,  making  the  column  behave  'softer'  than  predicted. 

Note  that  all  these  facts  have  direct  influence  on  the  column  stiffness,  explaining 
the  difference  between  the  FLPIER  model  and  the  test  results. 

A  'pseudo-dynamic'  test,  following  the  same  idea  presented  in  Example  1,  but 
with  an  adjustment  in  the  mass  and  damper  because  of  the  change  in  stiffness,  was  then 
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run  to  compare  the  effects  of  including  confinement  in  the  analysis.  A  plot  comparing 
the  FLPIER  unconfined  and  confined  dynamic  models  to  the  reference  is  shown  in  Fig. 
7-8.  Note  that  there  is  still  some  discrepancy  in  the  predicted  and  test  stiffness  of  the 
column,  which  may  be  caused  by  the  facts  mentioned  earlier.  Although  all  the 
approaches  could  predict  very  well  the  ultimate  moment,  only  the  confined  model  could 
predict  the  more  ductile  behavior  shown  in  Fig.  7-8  for  the  original  test. 

This  test  is  important  because  it  shows  clearly  the  effects  of  confinement.  Note 
that  for  the  static  analysis  the  theoretical  column  strength  is  just  slightly  over  the  actual 
test  value,  however  the  maximum  displacement  achieved  is  much  less  than  the  one  in 
the  test.  This  is  because  no  confinement  was  considered.  This  illustrates  the  main 
characteristic  of  confinement,  that  although  there  is  no  great  increase  in  the  column 
strength  (Fig.  7-8),  the  column's  ductility  is  reasonably  increased,  an  important 
characteristic  that  should  be  present  in  columns  in  seismic  regions.  Although  the 
confined  model  gives  a  response  much  closer  to  the  original  test,  failure  of  the  column 
at  about  23  in  can  not  be  predicted. 
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Fig.  7-6.  Shear  comparison  Example  3  x  FLPIER 
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Chai  et  ah- 
FLPIER  secant-|-- 
FLPIERtangeivQ- 


10  15 
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Fig.  7-7.  Moment  comparison  Example  3  and  FLPIER 
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Fig.  7-8.  Dynamic  model  comparison 
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Example  4  -  Circular  Reinforced  Concrete  Column  2 

Another  mono  tonic  test  for  a  circular  column  was  presented  in  Chai  et  al.  (1991). 
The  column  test  data  is  summarized  in  Table  7-4. 

Table  7-4.  Design  details  for  Example  4 


Diameter 

24" 

Height 

12' 

Cover  to  main  bar 

0.8" 

Concrete  strength  fc 
Modulus  of  elasticity  Ec 

4.725  ksi 
3918  ksi 

Longitudinal  steel 

Yield  strength^ 

Modulus  of  elasticity  Es  (adopted) 

26#6 
45.7  ksi 
29000  ksi 

Transverse  steel 
Yield  strength^ 

#2  hoop  at  5" 
51  ksi 

Axial  force 

400  kips 

Following  the  same  procedure  used  in  Example  3  an  incremental  load  static  test 
was  first  performed  using  the  secant  and  tangent  approaches.  No  confinement  was 
considered.  The  comparison  for  the  shear  forces  can  be  seen  in  Fig.  7-9.  Note  the  effect 
of  the  second  order  moment  as  the  shear  force  decreases  after  2  in  of  lateral 
displacement.  Then  a  pseudo-dynamic  test  was  run  to  verify  the  FLPIER  predictions,  for 
the  confined  and  unconfined  options.  The  comparison  can  be  seen  in  Fig.  7-10  and  Fig. 
7-11.  Note  that  we  now  have  a  much  closer  agreement  in  the  stiffness.  This  agreement 
can  be  explained  by  the  fact  that  this  is  a  much  smaller  column,  so  the  factors  that 
contributed  to  the  discrepancy  in  the  previous  test,  have  a  much  smaller  influence  in  this 
test. 

It  is  interesting  to  note  that  based  on  the  FLPIER  confinement  model  the  column 
is  not  adequately  confined.  Note  that  in  Fig.  7-10  the  confined  response  is  even  worse 
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than  the  unconfined  response.  The  conclusion  is  that  if  the  column  is  not  adequately 
confined  the  confinement  model  is  going  to  give  the  worse  results,  and  therefore  should 
not  be  used.  However  when  the  hoop  spacing  was  changed  to  2  in,  a  much  closer 
agreement  can  be  observed  between  the  confined  model  and  the  actual  test  (Fig.  7-1 1). 
The  adopted  confined  model  makes  no  distinction  between  hoop  or  spiral  confinement. 
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Fig.  7-9.  Static  shear  comparison  Example  4 
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Fig.  7-10.  Dynamic  moment  comparison  Example  4.  Hoop  spacing  =  5  in 
lOOOOi 
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Fig.  7-11.  Dynamic  moment  comparison  Example  4.  Hoop  spacing  =  2  in 
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Example  5  -  Rectangular  Reinforced  Concrete  Column 

In  order  to  verify  the  concrete  model  under  cyclic  loading  a  series  of  tests 
performed  by  Bousias  et  al.  (1995)  was  used  for  comparison.  This  test  consisted  of  a 
series  of  1 1  tests  performed  to  study  the  behavior  of  reinforced  concrete  columns 
subjected  to  various  types  of  cyclic  loading.  The  tests  referred  as  SO,  SI,  S2,  S3,  S4  and 
S10  in  the  original  work,  were  used  as  comparison  because  they  reflect  more  closely  the 
type  of  cyclic  loading  common  in  dynamic  analysis.  Tests  SO  to  S4  consist  of  a 
cantilever  column  subjected  to  imposed  displacements  or  forces  in  the  X  and  Y  direction 
with  a  constant  axial  load.  Test  S10  is  also  a  cantilever  column,  but  the  axial  force  (Z 
direction)  is  not  constant  during  the  test.  These  directions  are  indicated  in  Fig.  7-12. 
Table  7-5  has  a  summary  of  the  loading  parameters  and  concrete  strength  for  each  test . 

Table  7-5.  Loading  for  each  test 


Test 

/c(Mpa) 

Axial  (KN) 

Loading  path 

SO 

30.75 

300 

Imposed  X  displacement 

SI 

29.0 

212 

Imposed  X  displacement 
Imposed  Y  displacement 

S2 

31.1 

284 

Imposed  X  displacement 

Imposed  Y  displacement 

S3 

29.9 

310 

Imposed  X  displacement 
Imposed  Y  force 

S4 

27.7 

253 

Imposed  X  displacement 
Imposed  Y  force 

S10 

28.5 

Variable 

Imposed  X  force 
Imposed  Z  force 

The  details  about  the  specimens  used  in  the  test  can  be  found  in  Gutierrez  et  al. 
(1993).  The  specimens  had  a  0.25-m-square  cross  section  and  a  free  length  of  1.50m, 
and  were  built  in  as  a  cantilever  into  a  1-m-square,  0.5  m  thick,  heavily  reinforced 
foundation  base.  Longitudinal  reinforcement  consisted  of  eight  16-mm-diameter  bars, 
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uniformly  distributed  around  the  perimeter.  A  double  hoop  pattern  of  8-mm-diameter 
stirrups  at  a  70  mm  spacing  was  used  as  transversal  reinforcement.  The  steel  yielding 
stress  is  460  MPa.  No  reference  is  made  in  the  original  work  about  the  values  of  the 
modulus  of  elasticity  for  steel  or  concrete. 

For  the  FLPIER  model,  at  first  16  nonlinear  discrete  elements  were  used  to 
model  the  column,  but  they  proved  to  be  very  time  consuming  for  these  cyclic  tests.  The 
number  of  discrete  elements  was  then  reduced  to  four,  and  a  good  agreement  could  still 
be  obtained.  All  the  moments  from  FLPIER  were  computed  at  the  internal  node  of 
element  four  indicated  in  Fig.  7- 1 2,  to  avoid  the  addition  of  the  second  order  moments, 
and  allow  the  results  to  be  compared  to  the  original  test.  The  original  tests  were  also 
setup  in  a  way  that  the  second  order  moments  were  avoided.  All  the  moments  reported 
are  about  the  Y  axis  (Fig.  7-12).  These  are  basically  cyclic  tests  with  loading  modeled  as 
imposed  tip  displacements,  so  a  lumped  mass  m,  a  damper  c,  and  a  spring  k  were  added 
to  the  top  of  the  pile  to  minimize  the  dynamic  effects.  The  computer  model  and  the 
values  for  the  dynamic  parameters  can  be  seen  in  Fig.  7-12. 
Monotonic  Tests 

The  first  step  in  the  computer  analsyis  was  to  verify  the  concrete  model  under 
monotonic  load.  Using  the  actual  static  version  of  the  computer  program  FLPIER,  the 
column  capacity  was  found.  This  was  done  by  applying  incremental  increasing  forces 
until  no  convergence  was  achieved.  Then  with  the  dynamics  option  on,  the  same  type  of 
analysis  was  done.  The  plot  of  moment  at  the  base  of  the  pile  versus  tip  displacements 
for  both  analysis  is  illustrated  in  Fig.  7-13.  Note  that  the  static  version  can  not  predict 
the  softening  behavior  of  the  column  beyond  failure  because  the  secant  stiffness 
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approach  is  used.  In  the  dynamics  version  the  tangent  stiffness  approach  is  used,  so  the 
softening  behavior  can  be  modeled.  The  moment  capacity  for  this  column,  from  the 
static  analsyis,  was  found  to  be  approximately  100,000  KN.mm,  at  a  correpondent  tip 
displacement  of  15.4  mm.  Note  from  Fig.  7-13  that  both  models  are  in  very  good 
agreement. 


m  =  10  KN V/mm 


1=1500  mm 


k  =  1500  KN/mm 


c  =  245  KN s/mm 


Node  where 
moments  were 
computed 


93. 75  mm 


Fig.  7-12.  Computer  model  for  Test  3 
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Fig.  7-13.  Column  capacity  using  FLPIER 
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A  parametric  study  on  the  influence  of  changing  the  parameters  of  steel  and 
concrete  was  also  performed  for  the  monotonic  case.  Each  of  the  parameters,  modulus 
of  elasticity  of  concrete  Ec,  compressive  strength  of  concrete /c,  modulus  of  elasticity  of 
steel  Es,  and  yield  stress  of  steel  fy,  was  changed  by  a  factor  of  0.5  for  comparison. 
Table  7-6  helps  to  identify  the  tests  and  the  changed  parameters.  The  input  files  are 
identified  by  the  .in  extension.  The  extension  .DAT  identifies  the  output  data  file  for  the 
plots.  Confinement  was  not  considered  in  these  tests.  The  comparison  between  the 
various  parametric  tests  and  the  model  with  the  original  parameters  is  shown  in  Figs.  7- 
14  to  Fig.  7-17. 


Table  7-6.  Parametric  tests,  units  are  KN/mm2 

File 

Ec 

fe 

Es 

fy 

Confinement 

bml.in 

26.25 

0.031 

200 

0.46 

No 

bml_l.in 

13 

0.031 

200 

0.46 

No 

bml_2.in 

26.25 

0.031 

100 

0.46 

No 

bml_3.in 

26.25 

0.015 

200 

0.46 

No 

bml_4.in 

26.25 

0.031 

200 

0.23 

No 

bml_5.in 

13 

0.031 

100 

0.46 

No 

bml_6.in 

26.25 

0.015 

200 

0.23 

No 

bml_7.in 

13 

0.015 

100 

0.23 

No 
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Fig.  7-14.  Ec  and  Es  changed 
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Fig.7-16.  £and/changed 
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Fig.  7-17.  All  parameters  changed 
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The  following  comments  can  be  made  about  these  monotonic  tests: 

a)  Note  that  when  the  modulus  of  elasticity  is  changed  there  is  no  significant 
decrease  in  the  column  strength,  note  however  that  reducing  the  modulus  of  elasticity  of 
steel  has  a  more  significant  influence  than  reducing  the  modulus  of  elasticity  of  concrete 
(Fig.  7-14). 

b)  From  Fig.  7-15  note  that  by  reducing  by  half  the  concrete  strength  decreases 
the  column  capacity,  but  not  in  the  same  proportion.  Also  note  that  by  reducing  the  steel 
strength  by  half  decreases  the  column  strength  by  a  larger  amount  than  in  the  concrete 
case. 

c)  From  Fig.  7-16  note  that  reducing  both  values  of  the  modulus  of  elasticity  by 
half  makes  the  model  unstable,  which  is  expected.  However  decreasing  the  strengths  by 
half,  decreases  the  column  capacity  by  approximately  the  same  proportion,  as  expected, 
but  the  model  did  not  became  unstable. 

d)  Finally  by  decreasing  all  the  parameters  by  half,  the  column  strength  and 
stiffness  is  decreased  by  half  as  expected. 

The  file  bml_5.in  did  not  converge.  This  is  the  test  where  the  elastic  modulus  of 
steel  and  concrete  were  decreased  by  half,  but  the  original  values  for  stresses  were 
maintained.  Note  that  decreasing  drastically  one  of  the  parameters,  while  keeping  the 
other  constant,  may  cause  instabilities  to  the  model. 

In  order  to  illustrate  the  strain  rate  effect,  another  set  of  tests  was  done 
considering  confinement  at  different  strain  rates.  The  strain  rates  were  changed  form  10" 
/s  (very  slow)  to  1/s  (very  fast).  The  concrete  and  steel  properties  were  not  changed, 
and  are  those  of  file  bml.in  as  described  in  Table  7-6.  Table  7-7  helps  to  identify  the 
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strain  rate  tests.  The  comparison  between  the  unconfined  and  confined  models,  as  well 
as  the  effect  of  different  strain  rates  can  be  seen  in  Fig.  7-18. 

Note  the  difference  between  the  unconfined  and  the  confined  model  and  the 
increase  in  the  moment  capacity  as  the  strain  rate  is  increased.  The  additional  strength 
provided  by  the  increase  in  the  strain  rate  should  be  used  with  caution.  In  the  case  of  an 
earthquake,  for  example,  it  is  very  probable  that  the  structure  will  experiment  some 
cracking  and  yielding,  what  will  make  it  less  stiff,  increasing  the  period  and  leading  to 
lower  strain  rates.  Therefore  the  use  of  small  strain  rates  is  a  safer  lower  bound  solution. 

Table  7-7.  Parametric  tests  for  confinement  under  different  strain  rates 


File 

Strain  rate  (1/s) 

Confinement 

bml.in 

0 

No 

bml  srl.in 

1  x  10"5 

Yes 

bml  sr2.in 

1x10- 

Yes 

bml   sr3.in 

1  x  10"J 

Yes 

bml  sr4.in 

1  x  10"2 

Yes 

bml   sr5.in 

1  x  10" 

Yes 

bml  sr6.in 

1 

Yes 

10       20       30       40       50       60       70       80       90 
Top  lateral  displacement  (mm) 

Fig.  7-18.  Confinement  and  Strain  rate  effect 
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Cyclic  Test  SO 

In  order  to  verify  the  behavior  of  the  concrete  model  under  cyclic  loading, 
which  is  typical  of  dynamic  analysis,  a  cyclic  test  was  done  for  the  first  90  seconds  of 
the  original  test.  The  maximum  imposed  displacement  was  slightly  over  15mm,  which 
characterizes  failure  from  the  static  tests.  The  tip  displacement  history  for  the  first  90 
seconds  of  the  original  test  can  be  seen  in  Fig.  7-19.  It  is  also  compared  to  the  FLPIER 
history.  Note  that  they  are  in  perfect  agreement.  The  first  run  was  performed  with  the 
original  data  for  steel  and  concrete.  The  comparison  is  shown  in  Fig.  7-20.  Note  that  the 
FLPIER  model  is  stiffer  than  the  test.  It  was  then  necessary  to  adjust  the  values  of  the 
elastic  modulus  of  steel  and  concrete  to  calibrate  the  model.  In  lieu  of  information  from 
the  original  work,  the  initial  adopted  values  for  steel  and  concrete  were  respectively,  Es 
=  199.96  KN/mm2,  and  Ec  =  26.25  KN/mm2,  which  are  reasonable  values  for  steel  and 
concrete  (Meyer  (1996)).  However  a  good  agreement  between  the  test  and  FLPIER  was 
obtained  with  Es  =  90  KN/mm2,  and  Ec  =  20  KN/mm2.  Three  factors  may  have 
contributed  for  this  disagreement  between  the  test  and  FLPIER: 

a)  The  value  of  the  modulus  of  elasticity  for  concrete  was  computed  using  the 
empirical  formula  given  by  Meyer  (1996),  which  could  introduce  some  error. 

b)  The  modulus  of  elasticity  for  the  steel  was  adopted  considering  a  typical 
value  for  construction  steel.  This  could  also  introduce  some  error. 

c)  The  third,  and  maybe  the  most  important,  component  for  the  error  could  be 
anchorage  slip.  Alsiwat  and  Saatcioglu  (1992),  and  Saatcioglu  et  al.(1992)  reported  that 
anchorage  slip  is  one  of  the  major  components  of  inelastic  deformation  in  reinforced 
concrete.  Alsiwat  and  Saatcioglu  (1992,  pp.  1)  define  anchorage  slip  as: 
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Anchorage  slips  occurs  when  the  critical  section  of  a  member  for  flexure  is 
located  near  the  adjoining  member.  Formation  of  a  flexural  crack  at  the  interface 
of  two  members  strains  the  reinforcement  crossing  the  crack.  Widening  of  the 
crack  may  produce  inelastic  strains  in  the  reinforcement.  This  results  in 
penetration  of  yielding  into  the  adjoining  member,  giving  rise  to  significant 
extension  of  reinforcement.  Additional  rigid  body  deformation  may  occur  due  to 
slippage  of  reinforcement.  The  combined  effect  of  reinforcement  extension  and 
slip  in  the  adjoining  member  may  be  referred  to  as  anchorage  slip. 

Anchorage  slip  will  result  in  member  end  rotations  that  are  not  accounted  for  in 

flexural  analysis,  resulting  in  a  member  that  is  'softer'  than  initially  predicted.  In  the 

FLPIER  proposed  concrete  model  this  effect  is  not  considered,  but  it  seems,  based  on 

the  results  presented  in  this  work,  that  this  effect  may  be  considered  by  decreasing  the 

modulus  of  elasticity  of  the  materials,  making  the  member  'softer'.   Alsiwat  and 

Saatcioglu  (1992),  and  Saatcioglu  et  al.  (1992)  developed  models  for  anchorage  slip, 

and  the  reader  is  referred  to  their  work  for  more  information  on  this  topic.  The  strain 

rate  adopted  for  all  the  confined  tests  was  10"05l/sec  (very  slow). 


Fig.  7-19.  Imposed  displacement  history  for  the  first  90  seconds 
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The  next  step  was  to  run  the  complete  test  changing  the  modulus  of  elasticity  for 
steel  and  concrete.  The  total  duration  of  the  test  is  700  seconds.  The  complete  imposed 
tip  displacement  history  is  shown  in  Fig.  7-22.  It  is  important  to  note  that  failure  for  the 
column  is  characterized  at  about  15  mm.  Note  that  in  the  test  the  column  is  subjected  to 
tip  displacements  of  about  90  mm,  what  is  very  extreme.  The  analysis  was  done  with  a 
time  step  equal  to  0.01s.  Table  7-8  helps  to  identify  the  modulus  of  elasticity  of  the 
materials  for  each  run.  The  confinement  was  also  added  to  the  analysis.  Figures  7-23  to 
7-26  show  the  comparison  of  each  model  to  the  test  results.  Note  that  the  response  given 
by  FLPIER  using  the  original  material  properties  is  always  stiffer  than  the  original 
response,  although  the  column  strength  is  very  close  for  all  cases.  However  when  the 
modulus  of  elasticity  for  steel  and  concrete  were  decreased,  note  that  there  is  a  much 
better  agreement,  but  the  FLPIER  model  is  still  stiffer.  It  seems  that  there  is  an 
additional  stiffness  degradation  factor,  caused  by  the  large  amplitude  of  the  cyclic 
loading,  not  accounted  for  by  the  analysis.  Note  that  cyclic  degradation  is  usually  related 
to  a  large  number  of  low  intensity  loading  cycles,  but  in  this  example,  after  a  relatively 
small  number  of  cycles,  the  stiffness  degradation  is  considerable.  This  may  be 
introduced  in  the  model  in  future  works. 


Table  7-8.  Parameters  for  test  SO-  units  are  KN/mm 


File 

tsOl 

ts02 

ts03 

ts04 

Ec 

26.25 

26.25 

20 

20 

Es 

200 

200 

90 

90 

Confinement 

No 

Yes 

No 

Yes 
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Fig.  7-22.  Imposed  tip  displacement  history  for  test 
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Fig.  7-25.  Comparison  ts03 
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Fig.  7-26.  Comparison  ts04 


Cyclic  Test  SI 

In  this  test  displacement  cycles  in  pairs  of  linearly  increasing  amplitude  are 
alternately  applied  in  the  two  transverse  directions  X  and  Y.  The  displacement  history 
for  each  direction  can  be  seen  in  Fig.  7-27  and  Fig.  7-28. 
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Fig.  7-27.  Imposed  displacement  in  X  direction 
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Fig.7-28.  Imposed  displacement  in  Y  direction 

The  axial  load  is  constant  and  is  equal  to  212  KN.  Again  a  parametric  study  was 
done  to  calibrate  the  model  until  the  best  possible  agreement  was  obtained.  Table  7-9 
has  all  the  parameters  for  the  tests.  The  values  for  the  concrete  and  steel  strength  are 
respectively,  fc  =  0.029  KN/mm2,  and  fy  =  0.46  KN/mm2.  The  comparison  between 
FLPIER  and  test  SI  is  shown  in  Fig.  7-29  to  Fig.  7-32.  The  moment  plotted  is  always 
about  the  y  axis,  and  the  top  lateral  displacement  is  always  in  the  X  direction,  according 
to  Fig.  7-12. 


Table  7-9.  Parameters  for  test  SI-  units  are  KN/mm2 


File 

tsll 

tsllc 

tsl2 

tsl2c 

Ec 

26.25 

26.25 

20 

20 

Es 

200 

200 

90 

90 

Confinement 

No 

Yes 

No 

Yes 
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Fig.  7-29.  Comparison  tsl  1 
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Fig.  7-32.  Comparison  tsl2c 
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Again  note  that  a  much  better  approximation  is  obtained  with  the  modified 
values  of  the  modulus  of  elasticity  of  the  materials.  Also  note  that  the  column  strength  is 
close  to  the  test  value  for  all  parameters,  however  note  that  the  more  ductile  column 
behavior,  without  such  significant  loss  in  strength,  is  only  obtained  with  the  confined 
model.  Note  that  in  Fig.  7-32  we  almost  have  a  perfect  match. 
Cyclic  Test  S2 

In  this  test  the  deflection  in  the  transverse  direction  X  is  held  constant,  while 
displacement-controlled  cycles  are  applied  in  the  orthogonal  direction  Y.  Note  that 
displacements  are  imposed  in  both  direction  so  another  spring  was  added  to  the  Y 
direction  at  the  top  of  the  column.  The  imposed  displacement  histories  in  the  directions 
X  and  Y  respectively  are  shown  in  Fig.  7-33  and  Fig.  7-34.  The  axial  load  is  284  KN. 
The  total  test  time  was  about  650  seconds,  the  adopted  time  step  was  0.22  seconds.  The 
material  data  is  fc  =  0.03 1  KN/mm2,  Ec  =  26.25  KN/mm2,/,  =  0.46KN/mm2,  and  Es  = 
200KN/mm2.  Confinement  was  not  used.  A  good  agreement  with  these  properties  was 
obtained.  The  comparison  can  be  seen  in  Fig.  7-35. 
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Fig.  7-33.  Imposed  displacements  in  X  direction 
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Fig.7-35.  Comparison  test  S2  x  FLPIER 
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Note  that  in  this  test  it  was  not  necessary  to  change  the  material's  modulus  of 
elasticity  for  a  good  agreement,  however  FLPIER  could  not  model  the  stiffness 
degradation  present  in  the  model  due  to  the  amplitude  of  the  imposed  displacements.  It 
is  interesting  to  note  that  the  same  stiffness  degradation  is  present  in  the  unloading.  Also 
note  that  the  column  strength  predicted  by  FLPIER  is  slightly  larger  than  the  one  given 
by  the  test.  This  indicates  that  like  the  stiffness  degradation  factor  mentioned  before,  a 
strength  degradation  factor,  based  on  the  displacement  amplitude,  rather  than  number  of 
cycles,  should  be  considered  in  the  model.  Observe,  however,  that  as  the  displacements 
get  larger,  the  moments  are  in  better  agreement,  what  indicates  that  the  softening  model 
works. 
Cyclic  Tests  S3  and  S4 

These  are  mixed  mode  control  tests:  displacement-controlled  in  the  X-direction 
and  force-controlled  in  the  Y-direction.  The  imposed  displacement  in  the  X-direction 
and  force  in  the  Y-direction  are  shown  in  Fig.7-36  and  Fig.7-37.  The  comparison 
between  the  test  and  FLPIER  is  shown  in  Fig.  7-38  and  7-39.  Test  S4  differs  form  S3  in 
that  each  constant  level  of  the  Y-force  is  applied  first  in  the  +Y  direction  and  then  in  the 
-Y  direction,  with  the  repetition  of  three  cycles  of  X-displacement.  The  imposed 
displacement  in  the  X-direction  and  force  in  the  Y-direction  for  test  S4  are  shown  in 
Fig. 7-41  and  Fig.7-42,  respectively.  The  comparison  is  shown  in  Fig.  7-43  and  7-44. 
The  test  parameters  are  given  in  table  7-10. 
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Table  7-10.  Parameters  used  in  test  S3  and  S4 


File 

Test 

fc 

Ec 

fy 

Es 

Confinement 

Time  step 

ts31.in 

S3 

0.029 

26.25 

0.46 

200 

No 

0.35 

ts32.in 

S3 

0.029 

20 

0.46 

90 

No 

0.35 

ts41.in 

S4 

0.027 

26.25 

0.46 

200 

No 

0.41 

ts42.in 

S4 

0.027 

20 

0.46 

90 

No 

0.41 

ts43.in 

S4 

0.027 

20 

0.25 

90 

No 

0.41 
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Fig.7-36.  Imposed  X  displacement  for  test  S3 
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Fig.7-37.  Imposed  Y  forces  for  test  S3 
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Fig.7-38.  Comparison  Test  S3  x  FLPIER  ts31 
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Fig.7-39.  Comparison  Test  S3  x  FLPIER  ts32 
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Fig.  7-40.  Imposed  displacement  in  X  direction  for  test  S4 
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Fig.  7-4 1 .  Imposed  load  in  Y  direction  for  test  S4 
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Fig.  7-44.  Comparison  ts43 


Note  again  that  by  changing  the  values  of  the  modulus  of  the  elasticity  of 
concrete  and  steel  improves  the  FLPIER  model,  however  FLPIER  again  over  predicted 
the  column  capacity,  which  also  indicates  the  necessity  of  a  strength  degradation  factor 
for  the  FLPIER  model.  In  ts43  the  steel  strength  was  reduced  just  to  illustrate  the  effect 
of  changing  one  of  material  parameters  in  the  cyclic  analysis. 
Cyclic  Test  S10 

Finally  S10  is  a  test  where  the  axial  load  (Z  direction)  is  variable.  Again  the 
computer  model  is  just  like  the  previous  examples,  with  a  mass  and  a  damper  attached  at 
the  top  of  the  column  to  minimize  the  dynamic  effects.  The  material  properties  used  for 
the  various  models  are  shown  in  Table  7. 1 1 .  The  imposed  force  in  the  X  direction  is 
shown  in  Fig.  7-45,  and  the  imposed  force  in  the  Z  direction  is  shown  in  Fig.  7-46. 
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The  results  for  the  various  tests  run  is  shown  in  Figs.  7-47  to  7-50. 


Table  7-11. 

Material  properties  for 

example  S10 

File 

A 

Ec 

fy 

Es 

Conf. 

A/(s) 

tslOl 

0.028 

26.25 

0.46 

200 

No 

0.1 

tsl02 

0.028 

20 

0.25 

90 

No 

0.1 

tsl03 

0.028 

26.25 

0.35 

200 

Yes 

0.1 

tsl04 

0.028 

20 

—2 

0.35 

90 

Yes 

0.1 
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Fig.  7-45.  Imposed  force  in  X-direction  for  S 10 
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Fig.  7-48.  Comparison  tsl02  x  test 
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Note  again  that  a  better  response  is  obtained  with  the  corrected  values  for  the 
modulus  of  elasticity  of  the  materials.  In  this  test  the  steel  strength  was  also  changed  to 
illustrate  its  effect  in  the  cyclic  response  of  the  column.  Note  from  Fig.  7-48  that  the 
column  strength  is  greatly  reduced  when  the  steel  strength  is  reduced.  In  this  test  the 
best  fit  was  obtained  using  a  steel  strength^  =  0.35  KN/mm2  and  the  confined  model, 
as  seen  in  Fig.  7-50.  This  discrepancy  in  the  steel  strength  may  indicate  that  the  steel 
used  for  this  column  may  be  weaker  than  the  specified  value.  It  may  also  indicate  that 
varying  the  axial  force  may  also  add  some  strength  degradation  to  the  column,  not 
accounted  for  by  the  actual  FLPIER  model.  The  data  for  the  original  set  seems  to  have  a 
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cut  off  for  the  maximum  and  minimum  values  for  the  moments  (flat  portion  of  test  data 
in  the  Figures),  but  the  original  work  does  not  mention  anything  about  it. 

The  following  paragraphs  summarize  the  observations  made  for  the  examples 
presented  in  this  section: 

a)  In  the  beginning  of  the  modeling  phase,  when  16  elements  were  included  in 
the  analysis,  the  system  became  unstable,  because  additional  hinges  formed  within  the 
additional  elements.  It  was  also  observed  that  the  model  is  very  sensitive  to  the  size  of 
the  time  step  used  in  the  analysis.  For  cycles  of  very  short  duration,  like  in  example  SO, 
very  small  time  steps  should  be  used  in  order  to  obtain  a  correct  response. 

b)  Although  reducing  the  values  of  the  modulus  of  elasticity  for  steel  and 
concrete  improves  the  model  behavior,  it  may  also  become  unstable  for  large 
displacements. 

c)  It  is  clear  that  for  larger  displacement  amplitudes  the  bilinear  model  is  not  a 
good  representation  of  the  steel  behavior.  Note  the  change  in  the  shape  of  the  hysteresis 
graphics  for  the  larger  amplitudes.  This  suggests  that  a  more  accurate  model  for  the 
steel,  including  strain  hardening  and  possibly  cyclic  degradation,  will  result  in  a  better 
response.  The  introduction  of  a  concrete  gap  seems  to  model  the  member  stiffness 
degradation  very  well. 

d)  The  changes  in  the  tangent  elastic  modulus  suggested  by  Soroushian  et  al. 
(1986),  were  not  used  since  it  is  was  found  that  the  model  is  very  sensitive  to  such 
changes. 

e)  The  gap  model  suggested  in  Chapter  3  was  abandoned  because  it  introduced 
instability  in  the  analysis. 
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f)  The  bilinear  model  was  not  tested  because  the  approximation  obtained  with 
the  rational  model  was  very  reasonable. 

g)  It  is  essential  that  the  effects  of  anchorage  slip  be  included  in  the  nonlinear 
cyclic  analysis  of  reinforced  concrete  members  subjected  to  extreme  cyclic  loading. 

Example  6  -  Piles  in  Sand 

In  this  example  the  response  given  by  FLPIER  considering  the  soil  effects  is 
compared  to  tests  conducted  by  Wilson  et  al.  (1997).  This  was  a  series  of  dynamic 
centrifuge  tests  to  investigate  soil-pile-structure  interaction  in  liquefiable  sand.  The 
models  consisted  of  structures  supported  on  single  piles,  and  2x2  and  3x3  pile  groups. 
Although  in  the  original  work  a  variety  of  tests  were  performed,  three  tests  were  chosen 
for  comparison.  From  the  series  of  tests  referred  as  CSP2,  a  single  pile  with  a  single 
column  called  SP,  a  2  x  2  pile  group  with  single  a  column  called  PG2,  and  a  3  x  3  pile 
group  with  single  column  called  PG3,  were  selected  for  comparison.  In  the  original  test 
the  models  were  subjected  to  a  series  of  shaking  events,  with  liquefaction  occurring 
during  the  stronger  events.  In  this  work,  to  avoid  the  effects  of  liquefaction  not  yet 
included  in  the  model,  the  structure  is  subjected  only  to  the  first  event  on  the  series,  and 
the  results  are  then  compared.  The  structure  is  assumed  to  behave  linearly  in  all  the 
tests.  The  three  tests  are  described  next. 
Test  SP 

System  SP  simulates  a  superstructure  with  a  mass  of  50  tonnes  (1  tonne  =  1000 
kg)  supported  by  a  single  steel  pipe  pile  0.67  m  in  diameter,  16.8  m  long,  and  with  a  19 
mm  wall  thickness.  The  column  height  is  5.4  m.  Figure  7-51  shows  the  model  used  in 
FLPIER.  The  soil  data  for  all  the  tests  is  described  in  Table  7-12.  A  set  of  five  tests  were 
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run,  where  only  the  shear  wave  velocity  for  the  soil  layers  was  changed.  Because  the 
values  for  the  shear  wave  velocity  for  the  soil  were  not  given  in  the  original  work, 
reasonable  values  obtained  from  Richart  et  al.  (1970)  were  adopted.  The  values  vi  =  68 
m/s  and  \2  =  78  m/s  were  adopted  for  the  upper  and  lower  layers  respectively,  because 
they  gave  the  best  response.  The  values  of  the  soil  dampers  are  based  on  these  values  as 
mentioned  on  Chapter  6.  No  damping  was  added  to  the  structure.  The  acceleration 
record  for  the  top  ring  of  the  centrifuge  was  taken  as  input  for  the  FLPIER,  it  can  be 
seen  in  Fig.  7-52.  The  comparison  between  FLPIER  and  the  test  for  displacements  and 
moments  can  be  seen  in  Figs.  7-53  and  7-54  respectively.  Note  that  the  comparison  is 
very  good  for  a  single  pile. 


Table  7-12.  Some  soil  properties  for  CSP1 


Soil  Type 

Nevada  Sand 

pmax 

1.76 

pmin 

1.41 

Unit  Weight  (g/cmJ) 

1.52  Upper  layer 

1 .66  Lower  layer 

Relative  Density 

35-45  %  Upper  layer 

75-80%  Lower  layer 

After  Wilson  etal.  (1997) 

-j—  ]     M  =  50  tonnes 


5.4  m 


16.8  m 


Loose  saturated  Nevada  sand,  t  =  9. 1  m 


Dense  saturated  Nevada  sand,  t  =  7.6  m 
Fig.  7-51.  Single  pile  in  sand 
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Fig.  7-52.  Top  ring  acceleration 
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Fig.  7-53.  Top  displacement  comparison 
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Fig.  7-54.  Moments  comparison 


Tests  PG2  and  PG3 

Tests  PG2  and  PG3  consist  of  a  single  column  supported  by  a  2  x  2  (PG2)  and  a 
3x3  (PG3)  pile  group  respectively.  The  pile  and  soil  properties  are  identical  to  those  in 
the  single  pile  test  described  previously.  The  structure  dimensions  (in  meters)  and  setup 
are  shown  in  Figs.  7-55  and  7-56.  The  masses  for  the  cap,  column  and  top  mass  are 
given  in  Table  7-13.  The  period  of  the  superstructure  considering  it  is  fixed  at  the  base 
of  the  column  is  0.5  s.  The  prototype  model  in  FLPIER  had  its  pier  stiffness  adjusted  to 
have  this  same  period.  Because  the  structure  is  considered  to  be  linear  this  was  done 
simply  by  changing  the  values  of  the  modulus  of  elasticity  E,  and  moment  of  inertia  /, 
for  the  linear  3D  beam  elements  that  model  the  pier. 


Table  7-13.  Masses  (tons)  for  pile  group  tests 


PG 

Cap  mass 

Column  mass 

Top  mass 

2x2 

132 

174 

233 

3x3 

329 

193 

468 
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Fig.  7-55.  2x2  pile  group  in  sand 
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Fig.  7-56.  3x3  pile  group  in  sand 
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7-57.  Top  of  pile  lateral  displacement  comparison,  2x2  group 
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Fig.  7-58.  Top  of  pile  bending  moment  comparison,  2x2  group 
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Fig.  7-59.  Top  of  pile  bending  moment  comparison,  3x3  group 


The  comparison  between  FLPIER  and  the  test  is  shown  in  Figs.  7-57  and  7-58 
for  the  2  x  2  group  and  in  Fig.  7-59  for  the  3  x  3  group  (the  displacement  data  was  not 
available  for  this  configuration).  For  the  2  x  2  group  FLPIER  under  predicted  the 
displacements  as  it  can  be  seen  in  Fig.  7-57.  However  when  the  magnitude  of  these 
displacements  is  compared  to  the  structure  dimensions,  we  notice  that  they  are  very 
small.  Considering  that  we  are  trying  to  predict  a  maximum  displacement  of  about  0.002 
m  for  a  structure  that  is  about  30  m  (from  the  bottom  of  the  piles  to  the  top  of  the  piers), 
the  approximations  are  very  good.  Some  error  could  also  have  been  introduced  by 
incorrect  measurement  of  this  data. 
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On  the  other  hand  a  very  good  approximation  for  the  moment  was  obtained  from 
FLPIER  as  can  be  seen  in  Fig.  7-58.  For  the  3  x  3  group  such  a  good  approximation  was 
not  possible.  FLPIER  over  estimated  the  moments  by  about  50  percent.  Note  that  as  the 
number  of  piles  is  increased  the  results  for  the  moments  get  worse.  This  is  expected, 
since  the  group  influence,  not  included  in  the  FLPIER  model,  is  more  predominant.  It  is 
clear  that  more  study  is  necessary  in  this  area,  but  note  that  the  values  obtained  for  the 
FLPIER  approximation  are  perfectly  good  for  a  pre-design.  The  advantage  of  using 
FLPIER  is  that  the  pile  and  structure  configurations  can  be  easily  changed,  and  once  the 
displacements  are  in  the  desirable  range  a  more  sophisticated  analysis  can  be  done.  In 
the  future,  with  a  larger  database,  the  FLPIER  soil  model  can  be  calibrated  to  give  even 
closer  approximations  for  earthquake  analysis. 

Example  7  -  Mississippi  Dynamic  Test 

The  last  test  is  one  of  the  few  real  scale  tests  done  to  analyze  the  dynamic 
response  of  a  pile  group.  The  actual  data  for  the  test  was  obtained  from  the  East 
Pascagula  River  test  program  report.  The  tests  details  can  be  found  in  Brown  (1998). 
Basically  the  tests  consists  of  six  prestressed  rectangular  piles  subjected  to  a  pulse  load. 
Four  of  the  piles  were  battered.  The  structure  is  shown  in  Figs.  7-60.  The  load  history, 
applied  to  the  piles  cap  as  illustrated  in  Fig.  7-60,  is  shown  in  Fig.  7-61.  The 
displacement  history  for  the  piles  cap  is  shown  in  Fig.  7-62.  The  test  results  were  then 
compared  to  the  FLPIER  program  using  different  sets  of  soil  properties.  The  first  and 
second  comparison  were  based  on  CPT  and  SPT  soil  properties  estimates.  In  the  third 
comparison  the  proposed  original  soil  properties  were  used. 
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Fig.  7-60.  Mississippi  test  structure 
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Fig.  7-61.  Load  history  for  Mississippi  test 
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Fig.  7-62.  Displacement  history  for  Mississippi  test 
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Fig.  7-63.  Comparison  test  x  FLPIER  with  CPT  soil  properties 


1.5 


1    0.5 


S    -0.5 


-1 


0.5 


2.5 


1  1.5  2 

time  (s) 
Fig.  7-64.  Comparison  test  and  FLPIER  with  SPT  soil  properties 
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Fig.  7-65.  Comparison  test  and  FLPIER  with  original  soil  properties 


From  the  comparisons  we  observe  that  the  FLPIER  models  are  in  good 
agreement  with  the  original  test  data,  specially  the  SPT  model.  Note  that  the  initial 
period  for  the  system  is  very  close  to  the  actual  test  in  all  the  cases,  what  shows  that  the 
estimate  for  the  initial  stiffness  of  the  system  is  very  good.  However  the  same  agreement 
is  not  observed  later  in  the  response.  Note  that  the  FLPIER  model  is  suffer  than  the 
actual  test.  This  is  expected  since  the  FLPIER  soil  model  does  not  include  soil 
degradation  yet.  Also  note  that  there  is  a  small  residual  displacement  at  the  end  of  the 
analysis,  which  shows  that  there  was  some  damage  to  the  piles.  This  is  also  found  in  the 
FLPIER  model,  showing  that  the  concrete  model  is  adequate  even  for  prestressed  piles. 
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It  is  expected  that  once  soil  degradation  is  added  to  the  analysis,  an  almost  perfect  match 
will  be  obtained.  It  is  interesting  to  note  that  FLPIER  gave  good  approximations  for  a 
rather  complicated  model  of  nonlinear  dynamic  analysis.  Note  the  presence  of 
prestressed  battered  piles  and  four  soil  layers.  However,  modeling  was  easy  and  straight 
forward  using  the  FLPIER  program. 


CHAPTER  8 
CONCLUSIONS 

A  model  for  the  nonlinear  dynamic  analysis  of  reinforced  concrete  members  has 
been  presented.  The  results  shown  in  Chapter  7  indicate  that  the  proposed  model  for 
concrete  is  in  reasonable  agreement  with  the  tests  results  reported.  The  fiber  modeling 
adopted  by  FLPIER  seems  to  be  effective  in  modeling  steel  as  well  as  circular  and 
square  reinforced  concrete  sections.  It  is  clear  from  the  test  results  that  anchorage  slip  is 
an  important  concern  when  analyzing  reinforced  concrete  members  under  cyclic 
loading,  typical  during  earthquakes.  However  it  seems  that  it  is  possible  to  model  the 
anchorage  slip  effects  by  decreasing  the  value  of  the  modulus  of  elasticity  for  concrete 
and  steel.  The  model  can  be  calibrated  by  performing  tests  for  different  cross  sections 
and  adopting  the  corrected  value  for  the  modulus  of  elasticity  in  the  analysis.  Note  that 
for  all  types  of  tests  performed  for  the  rectangular  cross  section  reported  in  Chapter  7, 
the  same  value  of  the  corrected  modulus  of  elasticity  for  steel  and  concrete  gave  better 
results,  showing  that  anchorage  slip  is  a  function  of  the  cross  section,  rather  than 
loading. 

It  is  important  to  notice  that  the  introduced  concrete  gap  seems  to  model  the 
stiffness  degradation  of  the  section  very  well.  However  it  is  clear  from  the  hysteresis 
charts  of  all  the  tests  that  the  bilinear  steel  model  is  not  adequate  for  larger  lateral 
displacements.  In  such  cases  a  steel  model  that  includes  strain  hardening  and  possibly 
cyclic  degradation  should  be  used.  Anchorage  slip  should  also  be  included  in  this  model. 
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In  the  case  of  the  structures  where  soil  was  present  it  is  clear  that  more  studies 
still  need  to  be  done  when  multiple  piles  are  used,  because  the  dynamic  group  effect  is 
not  yet  well  understood.  However  the  approximations  given  by  FLPIER  gave  reasonable 
estimates  for  the  displacements  and  forces  acting  in  the  structure  from  very  simple 
parameters. 

The  next  steps  in  the  research  are: 

1 )  Introduce  a  new  steel  constitutive  model  including  strain  hardening. 

2)  Study  the  behavior  of  full  structures  modeled  with  the  nonlinear  discrete 
element  and  compare  the  results  to  the  literature. 

3)  Study  the  behavior  of  other  cross  sections,  such  as  prestressed  sections, 
sections  with  voids  and  sections  with  steel  sections  or  pipes  encased. 

4)  Analyze  structures  modeled  with  the  nonlinear  discrete  element  including  the 
soil  effects  and  compare,  if  possible,  the  results  to  the  literature. 

5)  Perform  more  tests  with  pile  groups  to  calibrate  the  FLPIER  model  for  the 
pile  group  effects. 

The  approach  presented  will  require  more  study  to  be  accurate  for  more 
complicated  cases,  but  what  is  important  is  that  in  all  examples  the  program  gave  a  good 
idea  about  the  magnitude  of  the  forces  and  displacements  under  a  variety  of  extreme 
loads.  Note  for  example  the  case  of  the  pile  groups  in  sand.  At  first  it  seems  very 
difficult  to  predict  the  response  of  a  pier  under  earthquake  loads,  including  the  soil- 
structure  interaction,  nevertheless  FLPIER  gave  good  approximations  for  the  single  and 
2x2  pile  group,  and  for  the  3  x  3  group,  its  prediction  is  off  by  about  50%.  In  the  case 
of  the  Mississippi  example  note  that  FLPIER  gave  an  excellent  prediction  for  the 
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response,  even  with  complications  such  as  nonlinear  prestressed  battered  piles.  In  this 
manner  engineers  can  have  a  feeling  of  what  to  expect  from  a  structure  under  a  certain 
earthquake  or  dynamic  loading. 

The  nonlinear  dynamic  analysis  of  structures  is  a  very  complex  subjected,  and 
the  author  with  this  work  tried  to  show  a  fairly  simple  approximation  to  the  problem  that 
showed  good  results  for  the  types  of  structures  studied. 


APPENDIX  A 
MASS  UNITS 

Problems  in  structural  dynamics  are  based  on  Newton's  second  law,  that  is, 

Force  =  (Mass)(Acceleration)  Al 

and  these  quantities,  plus  other  directly  related  to  them,  must  be  expressed  in  a 
consistent  system  of  units.  At  the  present  time  engineering  practice  is  in  the  process  of 
conversion  from  English  engineering  units  (United  States  Customary  System)  to  the 
International  System  (SI).  In  English  engineering  units  dimensional  homogeneity  in  Eq. 
Al  is  obtained  when  the  force  is  given  in  Ibf,  the  mass  in  slugs,  and  the  acceleration  in 
ft/sec  .  This  is  the  English  ft-lbf-sec  system  of  units.  The  slug  is  a  derived  unit,  and  form 
Eq.Al: 

1  slug  =  1  lbf-sec2/ft  A2 

The  units  of  force  and  mass  usually  lead  to  confusion  because  of  the  use  of  the 
term  weight  as  a  quantity  to  mean  either  force  or  mass.  On  the  one  hand,  when  one 
speaks  of  an  object's  "weight"  ,  it  is  usually  the  mass,  that  is,  quantity  of  matter,  that  is 
referred  to.  On  the  other  hand,  in  scientific  and  technological  usage  term  "weight  of  a 
body"  has  usually  meant  the  force  which  applied  to  the  body,  would  give  an  acceleration 
equal  to  the  local  acceleration  of  free  fall.  This  is  the  "weight  that  would  be  measured  by 
a  spring  scale.  If  the  mass  of  a  body  is  given  in  pounds  (lbm)  it  must  be  divided  by  the 
acceleration  of  gravity  g  s  32.2  ft/sec2  to  obtain  the  mass  as  referred  to  in  Newton's 
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second  law,  that  is,  a  force  of  1  Ibf  is  exerted  on  a  mass  of  \lbm  by  the  gravitational  pull 
of  the  earth. 

An  other  source  of  confusion  relies  on  the  fact  that  usually  there  is  no  distinction 
between  the  Ibf  and  lbm,  they  are  simply  referred  as  lb.  Because  FLPIER  needs  the  mass 
density  of  the  elements  for  dynamic  analysis,  we  decided  to  include  the  procedure  of 
how  to  obtain  this  quantity  for  two  of  the  most  used  materials  in  construction:  concrete 
and  steel.  Assuming  first  that  in  the  English  system  g  is  given  by: 

g  =  32.2  ft/sec2  =  386  in/sec2  A3 

The  unit  weight  of  normal  weight  concrete  is  usually  taken  as  yc  =  150  Ibf/ft 
(note  the  units  of  Ibf),  therefore  the  mass  density,  p,  is  obtained  by  dividing  yc  by  the 
acceleration  of  gravity  g: 

p  =1^150  M£la4.7^--4.7^  A4 

g      32.2  ft3  ft  ft'  f 

Note  that  the  slug  is  defined  for  Ibf  and  ft.  If  the  units  are  kip  and  in  the 
following  should  be  done: 

p  =  i^^£l4  .  2.25*10-  «£  =  2.25x10-  &£  A5 

12"    ft*    in4  in"  in4 
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The  same  procedure  can  be  used  for  steel  with  ys  =  490  Ibf/ft3,  resulting  in: 


.  7.35x10" 


kip.s2 


A6 


The  same  confusion  happens  in  SI  units  because  the  old  unit  of  force,  the  kgf  (1 
kgf=  9.8  N),  is  still  used  to  designated  the  unit  weight  of  materials,  sometimes  without 
the  distinction  from  the  mass-kilogram  that  is  represented  by  kg.  Table  A-l  has  a 
summary  of  these  results  in  the  English  and  SI  units. 


Table  A-l.  Mass  density  units 


English 

Concrete 

Steel 

yObf/fi3) 

150 

490 

(•(.slug/ft3) 

4.67 

15.23 

p  (kip  s  /in4) 

2.25  xlO"07 

7.35  xlO"07 

SI 

Concrete 

Steel 

y(kgf/m3) 

2320 

7860 

p(KNs2/m4) 

2.32 

7.86 

APPENDIX  B 
GAUSS  QUADRATURE 

Quadrature  is  the  name  applied  to  evaluating  an  integral  numerically,  rather  than 
analytically  as  is  done  in  tables  of  integrals.  There  are  many  quadrature  rules,  such  as 
Newton-Cotes  or  Simpson's  rule.  Here  we  discuss  only  the  Gauss-Legendre  rules, 
because  it  is  the  most  used  in  Finite  Elements  analysis.  A  brief  explanation  of  the 
method  extracted  from  Mathews  (1987)  follows. 

We  wish  to  find  the  area  under  the  curve 

y  =  f(x),        -1<x<1.  Bl 

What  method  gives  the  best  answer  if  only  two  function  evaluations  are  to  be 
made?  The  trapezoidal  rule  is  a  method  for  finding  the  are  under  the  curve,  it  use  two 
function  evaluations  at  the  end-points  (-\fi-\)),  (\J{\)).  But  if  the  graph  y  =  fix)  is 
concave  down,  the  error  in  approximation  is  the  entire  region  that  lies  between  the  curve 
and  the  line  segment  joining  the  points  (see  Fig.  B-l(a)). 
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t     y=A*) 


y=A*) 


(a) 


(b) 


Fig.  B-l.  (a)  Trapezoidal  approximation  using  the  abscissas  -1  and  1.  (b)  Trapezoidal 
approximation  using  abscissas  x;  and  x2  -  after  Mathews,  1987 


If  we  can  use  nodes  xi  and  x2  that  lie  inside  the  interval,  the  line  through  the  two 
points  (x/,f(x;)),  (x2,J{x2))  crosses  the  curve  and  the  are  under  the  line  more  closely 
approximates  the  are  under  the  curve  (see  Fig.  B-l(b)).  The  equation  of  this  line  is 


?  =  /(*,)  + 


(*-x,)[/(x2  )-/(*,)] 


B2 


and  the  area  of  the  trapezoid  under  this  line  is 


area  =  -^-f(x1)--^--f(x2). 


B3 


Notice  that  the  trapezoidal  rule  is  a  special  case  of  Eq.  B2.  When  we  choose  x;  = 
-1,  X2  =  1,  and  h  =  2,  then 
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r(/,*)-|/(«,)~/(^)=*/(^)+/<*2)  •  B4 

We  shall  use  the  method  of  undetermined  coefficients  to  find  the  abscissas  xi,  %2 
and  the  weights  wi,  W2  so  that  the  formula 

[j{x)dx*wj(xt)  +  wj{x2)  B5 

is    exact    for    cubic    polynomials    (i.e.,     f(x)  =  a,xi+a1x2+alx  +  a0).    Since    four 

coefficients  w/,  wj,  xi,  and  X2  need  to  be  determined  in  Eq.  B5,  we  can  select  four 

conditions  to  be  satisfied,  using  the  fact  that  integration  is  additive,  it  will  suffice  that 

Eq.  B5  is  exact  for  the  four  functions  7fa)=l,  x,  x2,  x3.  The  four  integral  conditions  are 

i 
f{x)  =  \:\\dx  =  2  =  w,  +  w2.  B6 


B7 


f(x)  =  x  :  jxdx  =  0  =  w,x,  +  w2x2 . 

f(x)  =  x2  :  jx2dx  =  |  =  w,x,2  +  w2x2  .  B8 

-i 

i 
f(x)  =  Xs  :  jx3dx  =  0  =  w,x,3  +  w2x2  .  B9 

-i 

Solving  the  system  of  nonlinear  equations  indicated  above  we  obtain: 

w,  =  w2  =  1  and  -  x,  =  x2  =  J^  m  0.5773502692  .     BIO 

We  have  found  the  nodes  and  the  weights  that  make  up  the  two  point  Gauss- 
Legendre  rule.  Since  the  formula  is  exact  for  cubic  equations,  the  error  term  will  involve 
the  fourth  derivative. 
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The  general  Appoint  Gauss-Legendre  rule  is  exact  for  polynomials  of  degree  < 
2JV-1  and  the  numerical  integration  formula  is 

G(f,N)  =  wNJ(xNI)  +  w„2f(xN2)  +  ...  +  wNNf(xNK).  Bll 

The  abscissas  Xfi.k  and  weights  wyj+  to  be  used  have  been  tabulated  and  are 
easily  available;  Table  B-l  gives  the  values  up  to  eight  points.  Also  included  in  the  table 
is  the  form  of  the  error  term  E(fN)  that  corresponds  to  G(f,N),  and  it  can  be  used  to 
determine  the  accuracy  of  the  Gauss-Legendre  integration  formula. 

The  values  in  Table  B-l  in  general  have  no  easy  representation,  this  fact  makes 
the  method  less  attractive  for  human  beings  to  use  when  hand  calculations  are  required. 
But  once  the  values  are  stored  in  a  computer  it  is  easy  to  call  them  up  when  needed,  the 
nodes  are  actually  roots  of  the  Legendre  polynomial  and  the  corresponding  weights 
must  be  obtained  by  solving  a  system  of  equations.  The  general  formula  for  N  points 
can  then  be  written: 


lf(x)dx  =  fiwNkf(xNt)  +  E(f,N) 


B12 
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Table  B-l.  Gauss-Legendre  abscissas  and  weights 


N 

Abscissas  x/v* 

Weights  w,vi 

Truncation  error 
E(f,N) 

2 

-0.5773502692 
+0.5773502692 

1.000000000 
1.000000000 

135 

3 

±0.7745966692 
0.00000000000 

0.555555556 
0.888888888 

f(%) 
15,750 

4 

±0.8611363116 
±0.3399810436 

0.3478548451 
0.6521451549 

f{%) 
3,472,875 

5 

±0.9061798459 
±0.5384693101 

0.00000000000 

0.2369268851 
0.4786286705 
0.5688888888 

f0%) 

1,237,732,650 

6 

±0.9324695142 
±0.6612093865 
±0.2386191861 

0.1713244924 
0.3607615730 
0.4679139346 

/'2\c)2"[6\f 
[l  2!f  1 3! 

7 

±0.9491079123 
±0.7415311856 
±0.4058451514 
0.00000000000 

0.1294849662 
0.2797053915 
0.3818300505 
0.4179591837 

/(14>(c)2l5[7!]4 
[I4!fl5! 

8 

±0.9602898565 
±0.7966664774 
±0.5255324099 
±0.1834346425 

0.1012285363 
0.2223810345 
0.3137066459 
0.3626837834 

/(,6)(C)217[8!]' 
[l6(Pl7! 

Source:  Mathews  (1987). 

Theorem  (Interval  transformation  for  Gauss-Legendre  rules):  Let  xNik  and 

wNik  be  the  abscissa  and  weights  for  the  JV-point  Gauss-Legendre  rule  on  the  interval  -1 
<  x  <  1 .  To  apply  the  rule  to  integrate  fit)  ona<t<b,  use  the  change  of  variables 


a  +  b       b-a 
t  = +  x 

2  2 


and 


b-a  , 
dt  = dx. 

2 


B13 


Then  the  integral  relationship  is 


//(,)*-//  ^  +  x 


a+b       b-a\b-a 


2    J    2 


-dx 


B14 


and  the  Gauss-Legendre  rule  for  [a,b]  is 
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J/(0*  =  —J]  "V.J1 —  +  V*  —J  ■  B15 

Two  and  Three  Dimensions.  Multidimensional  Gauss  rules,  called  Gaussian 
product  rules,  are  formed  by  successive  application  of  one-dimensional  Gauss  rules.  In 
two  dimensions,  consider  the  functionyfXjy),  the  formula,  neglecting  the  error  term,  is 

I  [f(x,y)dxdy  =  YZ^NJwNkf(xNJ)f(xNk)  B16 

and  in  three  dimensions,  considering  the  function  flxj/j)  the  formula  becomes 

lii  N     n     n 

I, 1,  U^,y,z)dxdyd2  =  2EZH,«.'vvivJ^,*/(^,i)/Kj)/W.J)  •      B17 

i=i  y=i  *=i 

It  is  not  necessary  to  use  the  same  Gauss  rule  in  all  directions,  but  doing  so  is 
most  common. 


APPENDIX  C 
FLPIER  MANUAL  FOR  DYNAMIC  ANALSYIS 


The  Flpier  program  now  supports  two  types  of  stiffness,  Secant  stiffness  and  Tangent  stiffness,  meaning  the  way  the 
stiffness  of  a  point  in  the  section  is  obtained  in  the  stress-strain  diagram. 

STIF 
S=NSTIF 

where       NSTIF      is  the  type  of  stiffness. 

NSTIF  =   0  -  secant  stiffness  (default). 

NSTIF  =    1  -  tangent  stiffness.  This  option  must  be  used  for  nonlinear  steel  sections  with  a  time  step  analysis. 

For  static  analysis  it  may  converge  faster  than  the  secant  method. 

This  section  MUST  end  with  a  blank  line. 

Type  of  Analysis 

There  are  two  options  for  the  type  of  analysis.  Condensed  means  that  all  the  structure  will  be  condensed  to  the  top  of 
the  piles,  and  then  the  problem  will  be  solved.  Full  means  that  there  will  be  no  condensation  for  the  analysis.  Because 
the  mass  condensation  is  not  exact  like  the  stiffness  condensation,  use  the  full  option  for  dynamic  step  by  step  analysis 
when  there  is  a  structure.  For  single  pile  analysis  use  the  condensed  option. 

SCND 
B=NFULL 

where       NFULL      is  the  type  of  analysis. 

NFULL  =  0  -  condensed  analysis  (default). 
NFULL  ■  1  -  full  analysis. 

This  section  MUST  end  with  a  blank  line. 

Dynamics 

For  dynamic  analysis  the  flag  showed  below  must  be  present  in  the  INPUT  file.  The  next  line  must  have  the  information 
described  below 

DYN 

Y=NDYNS  C=NDAMP  F=ALPHA1,BETA1,ALPHA2,BETA2  S=NPMAX  J=DMP  K=DMS  ONPRT  M=SMASS 

H=NSHM  L=NBMM  N=NDYSOL  U=NMSE  R=NCMOD  P=D1,D2,D3,D4,D5,D6,D7 

where       NDYNS    is  the  type  of  dynamics  solution. 

NDYNS  =0  -  Step  by  step  integration  (default). 

NDYNS  =  1  -  Spectrum  analysis  (for  the  structure  only). 
For  time  step  analysis  you  have  to  input  all  data  below,  for  spectrum  analysis  you  have  to  input  only  the  mass 
density  for  the  structure,  the  rest  of  the  data  is  specified  under  the  SPECTRUM  label.  If  there  is  no  reference 
in  the  input  file  for  any  of  the  variables,  the  value  assigned  for  the  variable  is  zero. 
NDAMP    is  the  damping  option. 

NDAMP  =0  -  no  damping  (default). 

NDAMP  =1  -damping. 
ALPHA1  ,BETA1       coefficients  for  Rayleigh  damping  for  the  structure  (C=ALPHA*M+BETA'K). 
ALPHA2,BETA2       coefficients  for  Rayleigh  damping  for  the  piles. 
ALPHA3,BETA3       coefficients  for  Rayleigh  damping  for  the  soil. 
NPMAX    is  the  maximum  number  of  time  steps  for  the  analysis. 
DMP         is  the  mass  density  for  the  piles. 
DMS         is  the  mass  density  for  the  structure. 
NPRT       is  the  output  option  for  time  step  analysis. 
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NPRT  =    0  -  maximum  displacements  and  maximum  forces  caused  by  the  maximum 
displacements  (default). 

NPRT  =    1  -  all  displacements  and  maximum  forces. 

NPRT  =    2  -  maximum  displacements  and  all  forces. 

NPRT  =    3  -  all  displacements  and  all  forces. 

Note  that  NPRT  =  2  or  3  options  only  allow  the  program  to  compute  the  element  forces  for  the 
options  above  (because  this  may  take  some  time  for  a  large  structure),  to  print  them  out  you  still  have  to  use  these  in 
addition  to  the  print  out  option.  For  example,  if  you  want  the  pile  forces  for  every  time  step  use  0=4  and  set  P=1  under 
the  PRINT  label.  If  you  want  only  the  structure  forces  use  T=1 .  The  maximum  forces  are  the  forces  caused  by  the 
maximum  displacements,  note  that  these  can  be  smaller  than  the  maximum  forces  for  the  structure.  For  options  2  and  3 
a  summary  of  the  maximum  forces  and  the  time  step  when  it  occurred  will  be  printed  out  at  the  end. 

SMASS    is  the  concentrated  mass  adopted  for  the  soil.  This  mass  is  applied  to  all  the  translational  DOF  X,  Y 

and  Z  to  represent  the  attached  soil  mass. 
NSHM      is  the  option  for  the  mass  matrix  for  the  cap. 

NSHM  =  0  -  consistent  mass  matrix  (default). 

NSHM  =  1  -  lumped  mass  matrix. 
NBMM      is  the  option  for  the  mass  matrix  for  the  structure  and  piles. 

NBMM  =  0  -  consistent  mass  matrix  (default). 

NBMM  =  1  -  lumped  mass  matrix. 
NDYSOL  is  the  option  for  the  type  of  numerical  solution.  This  option  is  only  valid  for  step  by  step  solution 
NDYNS=0 

NDYSOL  =  0  -  Newmark's  method  -  average  acceleration  (default). 

NDYSOL  =  1  -  Newmark's  method  -  linear  acceleration  (default). 

NDYSOL  =  2  -  Wilson-Theta  method. 
NMSE       is  the  option  for  multiple  support  excitation. 

NMSE  =  0  -  standard  analysis. 

NMSE  =  1  -  multiple  support  excitation. 
NCMOD   is  the  option  for  the  concrete  model  in  nonlinear  analysis. 

NCMOD  =  0  -  rational  model  1 . 

NCMOD  =  1  -  rational  model  1  with  crushing  and  cracking  option  ON. 

NCMOD  ■  2  -  rational  model  2. 

NCMOD=  3  -  rational  model  2  with  crushing  and  cracking  option  ON. 

NCMOD=  4  -  bilinear  model. 
D1...D7    this  is  an  option  if  the  user  wants  to  save  a  specific  NODE  displacement  to  a  file  for  possible  later 

plotting  or  checking. 

D1  -  is  the  number  of  NODES  that  will  me  saved  (maximum  =  6). 

D2...D7  -  is  the  number  of  the  NODE  to  be  saved.  For  example: 

P=6,1,2,3,4,5,6 

D1  =  6  -  six  NODES  will  be  saved  to  the  file. 

1...6-  NODES  1  to  6  will  be  saved  to  the  file.  They  do  not  have  to  be  in  any  specific  order. 

The  displacements  will  be  saved  in  a  text  file  with  the  name:  "inputname'.DSn,  the  velocities  in 
'inputname'.VSn,  and  the  accelerations  in  'inputaname'.ASn,  where  n  is  the  file  number  for  the  chosen  node.  In  the 
example  above  if  the  input  file  is  called  test. in,  the  displacements  for  node  1  are  saved  in  test.DSI,  for  node  2  in 
test  DS2  and  so  on. 

For  the  cap,  add  in  the  same  line  where  you  input  the  cap  properties  M  =  RMASS. 
where  RMASS  is  the  mass  density  for  the  cap  material. 

A  typical  input  line  for  step  by  step  analysis  would  look  like: 

DYN 

Y=0  C=1  F=.05,.Q5,.u,.u  S=100  J=1.309E-09  K=1.309E-09  0=1 

A  typical  input  line  for  the  spectrum  analysis  would  look  like: 

DYN 

Y=1  K=1.309E-09 

A  typical  input  line  for  the  cap  would  look  like: 

CAP 

E=4400  U=0.2  T=48  M=2.267E-07 

This  section  MUST  end  with  a  blank  line. 
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Spectrum  Analysis 

The  results  of  an  EIGEN  solution  can  be  used  to  perform  a  spectrum  analysis.  This  procedure  uses  the  mode 
superposition  method  to  combine  the  individual  eigen  vectors  into  a  single  response,  based  on  the  excitation  given  by  a 
response  spectrum.  Response  spectrums  are  usually  given  for  earthquake  loading.  This  procedure  combines  the 
individual  modes  response  for  a  spectrum  acting  in  each  of  the  three  global  directions,  X,Y  and  Z.  The  modal  responses 
are  combined  using  the  Complete  Quadratic  Combination  (CQC)  procedure,  the  directions  are  combined  using  Square 
Root  Sum  of  the  Squares  (SRSS).  SPECT  allows  only  a  single  input  response  spectrum  with  different  scale  factors  for 
that  spectrum  in  each  of  the  three  directions  (X,  Y,  and  Z).  The  program  will  allow  either  an  input  spectrum  or  an  input 
acceleration  record  and  generate  the  spectrum  values  internally. 

SPECTRUM 

For  input  Spectrum  use  the  following  lines: 

S=SP  D=DX,DY,DZ  N=NV 

where     SP  is  the  number  of  spectrum  points  used  to  define  the  response  spectrum  curve.  The  points  are  given 

in  pairs  (Circular  frequency  (rad  /  s),  Value)  and  are  assumed  linear  between  values. 

DX  is  the  scale  factor  to  apply  to  the  input  spectrum  for  use  in  the  X  direction.  The  response  spectrum 

values  are  scaled  by  this  factor  when  used  for  the  X  direction. 

DY  is  the  scale  factor  to  apply  to  the  input  spectrum  for  use  in  the  Y  direction.  The  response  spectrum 

values  are  scaled  by  this  factor  when  used  for  the  Y  direction. 

DZ  is  the  scale  factor  to  apply  to  the  input  spectrum  for  use  in  the  Z  direction.  The  response  spectrum 

values  are  scaled  by  this  factor  when  used  for  the  Z  direction. 

NV  is  the  number  of  eigenvectors  to  use  for  the  responses.  The  default  is  the  total  number  of  vectors 

solved  for  from  EIGEN. 

SPECTRUM  APPLIED  DIRECTIONS  (Repeat  for  as  many  lines  an  necessary) 

A  spectrum  analysis  also  needs  damping  ratios  for  each  mode  used  in  the  analysis.  The  damping  ratio  $si$  is  the 
percentage  of  damping  for  the  mode  in  question.  These  values  must  be  specified  if  a  spectrum  analysis  is  to  be 
performed.  The  next  lines  specifies  the  damping  ratios  to  use  for  each  mode.  As  many  lines  as  required  can  be  used. 

NF.NL.NI  S=S1 

where     NF  is  the  first  mode  in  a  generation  sequence  for  which  the  damping  ratio  is  used. 

NL  is  the  last  mode  in  a  generation  sequence  for  which  the  damping  ratio  is  used. 

Nl  is  the  increment  for  generating  mode  for  which  the  damping  ratio  is  used.  Modes  between  NF  and 

NL  will  also  use  the  specified  ratio. 

NL  and  Nl  can  be  left  blank  if  no  generation  is  desired. 

S1  is  the  damping  ration  value  to  be  used  (as  a  decimal). 

SPECTRUM  DEFINITION  LINES  (Repeat  for  as  many  lines  an  necessary) 

The  next  lines  specify  the  spectrum  function  values.  The  lines  MUST  contain  FOUR  pairs  of  numbers  each.  The 
number  of  lines  is  dictated  by  the  number  of  points  used  to  specify  the  function  (SP).  The  points  do  NOT  need  to  be  at 
even  spacing. 

F1,A1    F2,A2  F3.A3  F4.A4 

Where       Fi.Ai  are  the  frequency  and  acceleration  values  for  the  point  being  specified. 

NOTE:  the  last  point  in  the  spectrum  has  to  be  0,0. 
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Dynamic  Step  by  Step  Integration 

Any  structure  may  be  analyzed  using  the  step  by  step  integration  method  of  dynamic  analysis.  This  method  uses  the 
Newmark  method  and  Raleigh  damping  to  solve  for  dynamic  response  resulting  from  time  varying  loading.  The  time 
varying  loading  can  be  applied  as  a  single  load  function  applied  at  many  nodes  or  different  time  functions  applied  at 
specific  DOF.  The  STEP  program  will  perform  this  analysis.  The  distributed  and  concentrated  loads  applied  to  the 
structure  will  be  applied  as  a  CONSTANT  load  throughout  the  analysis.  The  applied  load  function  can  be  a  force  or  a 
ground  acceleration. 

TRANSIENT 

T-T1   L=L1   P=P1   Q=Q1   G=GF  A=A0  B=A1 

where     T1  is  the  time  step  increment  for  integration,  default  is  0.01s. 

L1  is  the  number  of  time  varying  load  functions  to  be  specified.  The  functions  are  applied  at  the 

specified  nodes  and  DOF. 

P1  is  the  maximum  number  of  time  points  specified  for  any  load  function.  If  three  functions  are 

specified,  P1  is  the  maximum  number  of  points  used  to  specify  any  of  the  three  functions. 

Q1  is  the  flag  for  the  type  of  load  function  applied.  If  Q1=0  then  the  load  is  a  force.  If  Q1=1,  then  the 

load  applied  is  a  ground  acceleration. 

GF  is  the  gravity  factor  to  multiply  times  the  acceleration  or  load  record  input  below. 

A0,  A1       are  global  damping  factors,  implying  that  only  one  damping  factor  will  be  applied  to  structure,  piles 
and  soil. 

Note:  you  must  choose  one  damping  factor,  either  here  or  under  the  DYN  header,  only  for  the  soil. 
Both  can  not  be  chosen. 

LOADING  FUNCTION  DEFINITIONS  (L1  Sets  of  lines) 

If  you  choose  the  multiple  support  excitation  option,  you  have  to  create  one  set  of  loading  functions  for  each  pile  node. 
Including  the  first  line. 

The  next  lines  specify  the  load  function  and  its  point  of  application.  There  should  be  L1  sets  of  lines. 

FIRST  LINES  :  Node  and  DOF  application  specification. 

There  can  be  as  many  of  these  lines  as  required  to  specify  all  loaded  nodes  and  DOF  for  this  load  function.  If  the  L= 
portion  is  NOT  specified,  ALL  active  DOF  will  be  loaded. 

If  Multiple  Support  Analysis  user  supplies  the  acceleration  record  for  each  pile  node,  1  to  1 6.  All  piles  are  assumed  to  be 
subjected  to  the  same  acceleration  record. 

N=NF,NL,NI  F=L1,L2,L3,L4,L51L6 

where     NF  is  the  first  node  in  a  generation  sequence  for  which  the  DOF  specification  is  used. 

NL  is  the  last  node  in  a  generation  sequence  for  which  the  DOF  specification  is  used. 

Nl  is  the  increment  for  generating  node  numbers  between  NF  and  NL  for  which  the  DOF  specification 

is  used.  NL  and  Nl  can  be  left  blank  if  no  generation  is  desired. 

For  Multiple  Support  Excitation  these  have  to  be  piles  nodes  limited  to  1  to  16. 

Li  is  the  state  at  which  the  ith  DOF  can  have,  either  loaded  or  NO  load.  Therefore  Li  can  have  ONLY 

the  two  following  values; 
Li  =  L  is  for  loaded. 
Li  =  N  is  for  NO  load. 

The  next  lines  specify  the  loading  function  values.  The  lines  MUST  contain  FOUR  pairs  of  numbers  each.  The  number 
of  lines  is  dictated  by  the  maximum  number  of  points  used  to  specify  the  function  (P1 ),  The  points  do  NOT  need  to  be  at 
even  spacing. 

T1.F1   T2.F2  T3.F3  T4.F4 


Where    Ti.Fi  are  the  time  and  force  (or  acceleration)  values  for  the  point  being  specified. 

This  section  MUST  end  with  a  blank  line. 

Point  Mass 

This  section  allows  the  addition  of  point  masses  to  a  structure. 

MASS 

The  next  line  specifies  the  mass  to  be  added  to  a  node  for  each  of  the  six  global  directions. 

NS.NF.NI  M=MX,MY,MZ,MRX,MRY,MR2 

where     NS  is  the  starting  node  to  add  the  mass  to. 

NF  is  the  final  node  to  add  the  mass  to. 

Nl  is  the  increment  to  generate  additional  node  numbers  at  between  NS  and  NF  at  which  to  add  mass. 

MX  \ 

MY  | 

MZ  }  are  the  mass  values  for  the  translational  and  rotational  X,Y,Z  directions 

MRX  | 

MRY  | 

MRZ  / 

This  section  must  end  with  a  blank  line. 

Point  Dampers 

This  section  allows  the  addition  of  point  dampers  to  a  structure. 

DAMP 

The  next  line  specifies  the  dampers  to  be  added  to  a  node  for  each  of  the  six  global  directions. 

NS.NF.NI  M=MX,MY,MZ,MRX,MRY1MRZ 

where     NS  is  the  starting  node  to  add  the  dampers  to. 

NF  is  the  final  node  to  add  the  dampers  to. 

Nl  is  the  increment  to  generate  additional  node  numbers  at  between  NS  and  NF  at  which  to  add 

dampers. 

MX  \ 

MY  | 

MZ  }  are  the  dampers  values  for  the  translational  and  rotational  X,Y,Z  directions 

MRX  | 
MRY 

MRZ  / 

You  can  not  add  concentrated  masses  or  dampers  to  the  pile's  nodes. 
This  section  must  end  with  a  blank  line. 
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Soil  Dynamic  Properties 

This  section  allows  the  user  to  specify  the  soil  properties  necessary  for  dynamic  analysis.  If  not  specified  they  are  set 
equal  to  zero. 

SDYN 

N=N50  T=SLR  V=SWVS1,SWVS2,...  F=SFDF 

where       N50  is  the  number  of  cycles  necessary  to  degrade  the  soil  by  50%. 

SLR  Is  the  rate  of  loading  for  slow  cyclic  loading. 

SWVSi  Is  the  shear  wave  velocity  for  each  soil  layer. 

SFDF  is  the  fully  degraded  soil  factor. 

This  section  must  end  with  a  blank  line. 

Free  Vibration  Option 

This  section  allows  the  user  to  run  a  free  vibration  problem  with  initial  imposed  parameters. 

FREEV 

T=DT  L=NUMN  Q=3  G-GF 

where      DT  is  the  time  step. 

NUMN  is  the  number  of  nodes  with  imposed  displacement  and  velocity. 

3  is  the  free  vibration  option 

GF  is  the  gravity  factor 

For  each  node,  there  must  be  a  line  with  the  following  information 

NODE  D  =  x  ,  y  ,  z  ,  9x  ,  8y  ,  0x  ,  xv  ,  yv,  zv,  Gxv  ,  Gyv  ,  6zv 

where        NODE  is  the  node  number  where  the  initial  conditions  are  applied. 

x.y.z  are  imposed  displacements  in  the  global  directions  X,Y,  and  Z  respectively. 

8x ,  By  ,  6x  are  imposed  rotations  about  the  global  axis  X,Y,  and  Z  respectively, 

xv  ,  yv,  zv  are  imposed  velocities  in  the  global  directions  X,Y,  and  Z  respectively. 

6xv  ,  6yv  ,  6zv  are  imposed  angular  velocities  about  the  global  axis  X ,  Y,  and  Z  respectively. 

Confinement 

This  section  allows  the  user  to  specify  the  parameters  for  the  concrete  confinement. 

DCONF 

F=  FYHOOP  S=  HOOPS  C=  HOOPC  D-  HOOPD  R=  EPP 

where       FYHOOP  is  the  yield  stress  for  the  hoop  steel. 

Note  for  English  units  use  KSI.  for  SI  units  use  MPA. 
HOOPS  is  the  hoop  spacing. 
HOOPC  is  the  hoop  core  diameter. 
HOOPD  is  the  hoop  bar  diameter. 
EPP  is  the  strain  rate  for  concrete. 
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